#This is the replication file for Thomann, E. van Engen, N. and L. Tummers. Forthcoming. The necessity of discretion: a behavioral evaluation of bottom-up implementation theory.  Journal of Public Administration Research and Theory.

# This code replicates the analysis of dataset 1, recoding method.

#R code based  on Thomann, Eva and Stefan Wittwer (2016). Performing fuzzy- and crisp set QCA with R: A user-oriented beginner's guide. Version 19.9.2016. URL: http://www.evathomann.com/links/qca-r-manual [date of access: 25.11.2016].


#clear workspace
rm(list = ls())


library(lattice); library(arm); library(xtable); library(foreign); library(psych); library(directlabels); library(betareg); library(VIM); library(base); library(plyr); library(dplyr); library(QCA); library(QCAGUI); library(SetMethods)

setwd("C:/Users/Eva/Dropbox/Arbeit/Publikationen/Papers/Two-factor theory/Material/Data/Dataset 1")


##############PREPARING THE DATASET###########

tfr <- read.csv("tfr_d1.csv", row.names=1, header = TRUE, sep = ";",  dec = ",")
describe(tfr)


#check missings

missings <- as.numeric(complete.cases(tfr))
tabmiss <- table(missings)
prop.table(tabmiss)

#identify  cases with missings on all indicators for a causal factor
##for W
colnames(tfr)

tfr$w1miss <- as.numeric(is.na(tfr$W1))
tabw1miss <- table(tfr$w1miss)
prop.table(tabw1miss)


tfr$w2miss <- as.numeric(is.na(tfr$W2))
tabw2miss <- table(tfr$w2miss)
prop.table(tabw2miss)

tfr$w3miss <- as.numeric(is.na(tfr$W3))
tabw3miss <- table(tfr$w3miss)
prop.table(tabw3miss)

tfr$w4miss <- as.numeric(is.na(tfr$W4))
tabw4miss <- table(tfr$w4miss)
prop.table(tabw4miss)

tfr$w5miss <- as.numeric(is.na(tfr$W5))
tabw5miss <- table(tfr$w5miss)
prop.table(tabw5miss)

wmiss <- subset(tfr, select = c("w1miss", "w2miss", "w3miss", "w4miss", "w5miss"))
tfr$wmiss <- rowSums(wmiss)
tfr$wmiss
tfr$wmissout <- as.numeric(tfr$wmiss ==5)
tfr$wmissout
tabwmissout <- table(tfr$wmissout)
prop.table(tabwmissout)


##for SP

tfr$sp1miss <- as.numeric(is.na(tfr$SP1))
tabsp1miss <- table(tfr$sp1miss)
prop.table(tabsp1miss)


tfr$sp2miss <- as.numeric(is.na(tfr$SP2))
tabsp2miss <- table(tfr$sp2miss)
prop.table(tabsp2miss)

tfr$sp3miss <- as.numeric(is.na(tfr$SP3))
tabsp3miss <- table(tfr$sp3miss)
prop.table(tabsp3miss)

tfr$sp4miss <- as.numeric(is.na(tfr$SP4))
tabsp4miss <- table(tfr$sp4miss)
prop.table(tabsp4miss)

tfr$sp5miss <- as.numeric(is.na(tfr$SP5))
tabsp5miss <- table(tfr$sp5miss)
prop.table(tabsp5miss)

tfr$sp6miss <- as.numeric(is.na(tfr$SP6))
tabsp6miss <- table(tfr$sp6miss)
prop.table(tabsp6miss)

spmiss <- subset(tfr, select = c("sp1miss", "sp2miss", "sp3miss", "sp4miss", "sp5miss", "sp6miss"))
tfr$spmiss <- rowSums(spmiss)
tfr$spmiss
tfr$spmissout <- as.numeric(tfr$spmiss==6)
tfr$spmissout
tabspmissout <- table(tfr$spmissout)
prop.table(tabspmissout)

##for TP
tfr$tp1miss <- as.numeric(is.na(tfr$TP1))
tabtp1miss <- table(tfr$tp1miss)
prop.table(tabtp1miss)


tfr$tp2miss <- as.numeric(is.na(tfr$TP2))
tabtp2miss <- table(tfr$tp2miss)
prop.table(tabtp2miss)

tfr$tp3miss <- as.numeric(is.na(tfr$TP3))
tabtp3miss <- table(tfr$tp3miss)
prop.table(tabtp3miss)

tfr$tp4miss <- as.numeric(is.na(tfr$TP4))
tabtp4miss <- table(tfr$tp4miss)
prop.table(tabtp4miss)

tfr$tp5miss <- as.numeric(is.na(tfr$TP5))
tabtp5miss <- table(tfr$tp5miss)
prop.table(tabtp5miss)

tfr$tp6miss <- as.numeric(is.na(tfr$TP6))
tabtp6miss <- table(tfr$tp6miss)
prop.table(tabtp6miss)


tpmiss <- subset(tfr, select = c("tp1miss", "tp2miss", "tp3miss", "tp4miss", "tp5miss", "tp6miss"))
tfr$tpmiss <- rowSums(tpmiss)

tfr$tpmissout <- as.numeric(tfr$tpmiss==6)
tfr$tpmissout
tabtpmissout <- table(tfr$tpmissout)
prop.table(tabtpmissout)

##for OP
tfr$op1miss <- as.numeric(is.na(tfr$OP1))
tabop1miss <- table(tfr$op1miss)
prop.table(tabop1miss)

tfr$op2miss <- as.numeric(is.na(tfr$OP2))
tabop2miss <- table(tfr$op2miss)
prop.table(tabop2miss)

tfr$op3miss <- as.numeric(is.na(tfr$OP3))
tabop3miss <- table(tfr$op3miss)
prop.table(tabop3miss)

tfr$op4miss <- as.numeric(is.na(tfr$OP4))
tabop4miss <- table(tfr$op4miss)
prop.table(tabop4miss)

tfr$op5miss <- as.numeric(is.na(tfr$OP5))
tabop5miss <- table(tfr$op5miss)
prop.table(tabop5miss)

tfr$op6miss <- as.numeric(is.na(tfr$OP6))
tabop6miss <- table(tfr$op6miss)
prop.table(tabop6miss)

opmiss <- subset(tfr, select = c("op1miss", "op2miss", "op3miss", "op4miss", "op5miss", "op6miss"))
tfr$opmiss <- rowSums(opmiss)

tfr$opmissout <- as.numeric(tfr$opmiss==6)
tfr$opmissout
tabopmissout <- table(tfr$opmissout)
prop.table(tabopmissout)

##for SM

tfr$sm1miss <- as.numeric(is.na(tfr$SM1))
tabsm1miss <- table(tfr$sm1miss)
prop.table(tabsm1miss)

tfr$sm2miss <- as.numeric(is.na(tfr$SM2))
tabsm2miss <- table(tfr$sm2miss)
prop.table(tabsm2miss)

tfr$sm3miss <- as.numeric(is.na(tfr$SM3))
tabsm3miss <- table(tfr$sm3miss)
prop.table(tabsm3miss)

tfr$sm4miss <- as.numeric(is.na(tfr$SM4))
tabsm4miss <- table(tfr$sm4miss)
prop.table(tabsm4miss)

tfr$sm5miss <- as.numeric(is.na(tfr$SM5))
tabsm5miss <- table(tfr$sm5miss)
prop.table(tabsm5miss)

tfr$sm6miss <- as.numeric(is.na(tfr$SM6))
tabsm6miss <- table(tfr$sm6miss)
prop.table(tabsm6miss)

tfr$sm7miss <- as.numeric(is.na(tfr$SM7))
tabsm7miss <- table(tfr$sm7miss)
prop.table(tabsm7miss)

tfr$sm8miss <- as.numeric(is.na(tfr$SM8))
tabsm8miss <- table(tfr$sm8miss)
prop.table(tabsm8miss)

tfr$sm9miss <- as.numeric(is.na(tfr$SM9))
tabsm9miss <- table(tfr$sm9miss)
prop.table(tabsm9miss)

tfr$sm10miss <- as.numeric(is.na(tfr$SM10))
tabsm10miss <- table(tfr$sm10miss)
prop.table(tabsm10miss)

tfr$sm11miss <- as.numeric(is.na(tfr$SM11))
tabsm11miss <- table(tfr$sm11miss)
prop.table(tabsm11miss)

tfr$sm12miss <- as.numeric(is.na(tfr$SM12))
tabsm12miss <- table(tfr$sm12miss)
prop.table(tabsm12miss)

smmiss <- subset(tfr, select = c("sm1miss", "sm2miss", "sm3miss", "sm4miss", "sm5miss", "sm6miss", "sm7miss", "sm8miss", "sm9miss", "sm10miss", "sm11miss", "sm12miss"))
tfr$smmiss <- rowSums(smmiss)

tfr$smmissout <- as.numeric(tfr$smmiss==12)
tfr$smmissout
tabsmmissout <- table(tfr$smmissout)
prop.table(tabsmmissout)

##for CM

tfr$cm1miss <- as.numeric(is.na(tfr$CM1))
tabcm1miss <- table(tfr$cm1miss)
prop.table(tabcm1miss)


tfr$cm2miss <- as.numeric(is.na(tfr$CM2))
tabcm2miss <- table(tfr$cm2miss)
prop.table(tabcm2miss)

tfr$cm3miss <- as.numeric(is.na(tfr$CM3))
tabcm3miss <- table(tfr$cm3miss)
prop.table(tabcm3miss)

tfr$cm4miss <- as.numeric(is.na(tfr$CM4))
tabcm4miss <- table(tfr$cm4miss)
prop.table(tabcm4miss)


cmmiss <- subset(tfr, select = c("cm1miss", "cm2miss", "cm3miss", "cm4miss"))
tfr$cmmiss <- rowSums(cmmiss)

tfr$cmmissout <- as.numeric(tfr$cmmiss ==4)
tfr$cmmissout
tabcmmissout <- table(tfr$cmmissout)
prop.table(tabcmmissout)

# calculate dropout rate due to missings

dropout <- subset(tfr, select = c("wmissout", "spmissout", "tpmissout", "opmissout", "smmissout", "cmmissout"))
tfr$dropout <- as.numeric(rowSums(dropout) == 6)
tabdropout <- table(tfr$dropout)
prop.table(tabdropout)

#delete cases with missings on all indicators for 1 causal factor

tfr = tfr[tfr$wmissout == 0,]
tfr = tfr[tfr$spmissout == 0,]
tfr = tfr[tfr$tpmissout == 0,]
tfr = tfr[tfr$opmissout == 0,]
tfr = tfr[tfr$smmissout == 0,]
tfr = tfr[tfr$cmmissout == 0,]
tfr

describe(tfr)

#check missings
missings <- as.numeric(complete.cases(tfr))
tabmiss <- table(missings)
prop.table(tabmiss)


#calculate dropout rate - from 1317 to 1004

do <- round(1-(1004/1317), digits=3)
do


#####CALIBRATION########

#calibrate indicator sets, using indirect calibration method
#the values on the likert scales for powerlessness and meaninglessness need to be reverted because they should mean powerfulness and meaningfulness
##indirect method

###for W
# rename raw variables
tfr <- rename(tfr, w1 = W1)
tfr <- rename(tfr, w2 = W2)
tfr <- rename(tfr, w3 = W3)
tfr <- rename(tfr, w4 = W4)
tfr <- rename(tfr, w5 = W5)
colnames (tfr)

tfr$W1 <- NA
tfr$W1 <- recode(tfr$w1, "5 = 1; 4 = 0.8; 3 = 0.4; 2 = 0.2; 1 = 0; else = copy")
describe(tfr$W1)

tfr$W2 <- NA
tfr$W2 <- recode(tfr$w2, "5 = 1; 4 = 0.8; 3 = 0.4; 2 = 0.2; 1 = 0; else = copy")
describe(tfr$W2)

tfr$W3 <- NA
tfr$W3 <- recode(tfr$w3, "5 = 1; 4 = 0.8; 3 = 0.4; 2 = 0.2; 1 = 0; else = copy")
describe(tfr$W3)

tfr$W4 <- NA
tfr$W4 <- recode(tfr$w4, "5 = 1; 4 = 0.8; 3 = 0.4; 2 = 0.2; 1 = 0; else = copy")
describe(tfr$W4)

tfr$W5 <- NA
tfr$W5 <- recode(tfr$w5, "5 = 1; 4 = 0.8; 3 = 0.4; 2 = 0.2; 1 = 0; else = copy")
describe(tfr$W5)

## for SP
# rename raw variables
tfr <- rename(tfr, sp1 = SP1)
tfr <- rename(tfr, sp2 = SP2)
tfr <- rename(tfr, sp3 = SP3)
tfr <- rename(tfr, sp4 = SP4)
tfr <- rename(tfr, sp5 = SP5)
tfr <- rename(tfr, sp6 = SP6)
colnames (tfr)

tfr$SP1 <- NA
tfr$SP1 <- recode(tfr$sp1, "5 = 0; 4 = 0.2; 3 = 0.4; 2 = 0.8; 1 = 1; else = copy")
describe(tfr$SP1)


tfr$SP2 <- NA
tfr$SP2 <- recode(tfr$sp2, "5 = 0; 4 = 0.2; 3 = 0.4; 2 = 0.8; 1 = 1; else = copy")
describe(tfr$SP2)


tfr$SP3 <- NA
tfr$SP3 <- recode(tfr$sp3, "5 = 0; 4 = 0.2; 3 = 0.4; 2 = 0.8; 1 = 1; else = copy")
describe(tfr$SP3)


tfr$SP4 <- NA
tfr$SP4 <- recode(tfr$sp4, "5 = 0; 4 = 0.2; 3 = 0.4; 2 = 0.8; 1 = 1; else = copy")
describe(tfr$SP4)


tfr$SP5 <- NA
tfr$SP5 <- recode(tfr$sp5, "5 = 0; 4 = 0.2; 3 = 0.4; 2 = 0.8; 1 = 1; else = copy")
describe(tfr$SP5)

tfr$SP6 <- NA
tfr$SP6 <- recode(tfr$sp6, "5 = 0; 4 = 0.2; 3 = 0.4; 2 = 0.8; 1 = 1; else = copy")
describe(tfr$SP6)

## for TP
tfr <- rename(tfr, tp1 = TP1)
tfr <- rename(tfr, tp2 = TP2)
tfr <- rename(tfr, tp3 = TP3)
tfr <- rename(tfr, tp4 = TP4)
tfr <- rename(tfr, tp5 = TP5)
tfr <- rename(tfr, tp6 = TP6)
colnames (tfr)

tfr$TP1 <- NA
tfr$TP1 <- recode(tfr$tp1, "5 = 0; 4 = 0.2; 3 = 0.4; 2 = 0.8; 1 = 1; else = copy")
describe(tfr$TP1)


tfr$TP2 <- NA
tfr$TP2 <- recode(tfr$tp2, "5 = 0; 4 = 0.2; 3 = 0.4; 2 = 0.8; 1 = 1; else = copy")
describe(tfr$TP2)


tfr$TP3 <- NA
tfr$TP3 <- recode(tfr$tp3, "5 = 0; 4 = 0.2; 3 = 0.4; 2 = 0.8; 1 = 1; else = copy")
describe(tfr$TP3)


tfr$TP4 <- NA
tfr$TP4 <- recode(tfr$tp4, "5 = 0; 4 = 0.2; 3 = 0.4; 2 = 0.8; 1 = 1; else = copy")
describe(tfr$TP4)


tfr$TP5 <- NA
tfr$TP5 <- recode(tfr$tp5, "5 = 0; 4 = 0.2; 3 = 0.4; 2 = 0.8; 1 = 1; else = copy")
describe(tfr$TP5)

tfr$TP6 <- NA
tfr$TP6 <- recode(tfr$tp6, "5 = 0; 4 = 0.2; 3 = 0.4; 2 = 0.8; 1 = 1; else = copy")
describe(tfr$TP6)


##for OP

tfr <- rename(tfr, op1 = OP1)
tfr <- rename(tfr, op2 = OP2)
tfr <- rename(tfr, op3 = OP3)
tfr <- rename(tfr, op4 = OP4)
tfr <- rename(tfr, op5 = OP5)
tfr <- rename(tfr, op6 = OP6)
colnames (tfr)

tfr$OP1 <- NA
tfr$OP1 <- recode(tfr$op1, "5 = 0; 4 = 0.2; 3 = 0.4; 2 = 0.8; 1 = 1; else = copy")
describe(tfr$OP1)


tfr$OP2 <- NA
tfr$OP2 <- recode(tfr$op2, "5 = 0; 4 = 0.2; 3 = 0.4; 2 = 0.8; 1 = 1; else = copy")
describe(tfr$OP2)


tfr$OP3 <- NA
tfr$OP3 <- recode(tfr$op3, "5 = 0; 4 = 0.2; 3 = 0.4; 2 = 0.8; 1 = 1; else = copy")
describe(tfr$OP3)


tfr$OP4 <- NA
tfr$OP4 <- recode(tfr$op4, "5 = 0; 4 = 0.2; 3 = 0.4; 2 = 0.8; 1 = 1; else = copy")
describe(tfr$OP4)


tfr$OP5 <- NA
tfr$OP5 <- recode(tfr$op5, "5 = 0; 4 = 0.2; 3 = 0.4; 2 = 0.8; 1 = 1; else = copy")
describe(tfr$OP5)

tfr$OP6 <- NA
tfr$OP6 <- recode(tfr$op6, "5 = 0; 4 = 0.2; 3 = 0.4; 2 = 0.8; 1 = 1; else = copy")
describe(tfr$OP6)


##for SM

tfr <- rename(tfr, sm1 = SM1)
tfr <- rename(tfr, sm2 = SM2)
tfr <- rename(tfr, sm3 = SM3)
tfr <- rename(tfr, sm4 = SM4)
tfr <- rename(tfr, sm5 = SM5)
tfr <- rename(tfr, sm6 = SM6)
tfr <- rename(tfr, sm7 = SM7)
tfr <- rename(tfr, sm8 = SM8)
tfr <- rename(tfr, sm9 = SM9)
tfr <- rename(tfr, sm10 = SM10)
tfr <- rename(tfr, sm11 = SM11)
tfr <- rename(tfr, sm12 = SM12)
colnames (tfr)

tfr$SM1 <- NA
tfr$SM1 <- recode(tfr$sm1, "5 = 0; 4 = 0.2; 3 = 0.4; 2 = 0.8; 1 = 1; else = copy")
describe(tfr$SM1)


tfr$SM2 <- NA
tfr$SM2 <- recode(tfr$sm2, "5 = 0; 4 = 0.2; 3 = 0.4; 2 = 0.8; 1 = 1; else = copy")
describe(tfr$SM2)


tfr$SM3 <- NA
tfr$SM3 <- recode(tfr$sm3, "5 = 0; 4 = 0.2; 3 = 0.4; 2 = 0.8; 1 = 1; else = copy")
describe(tfr$SM3)


tfr$SM4 <- NA
tfr$SM4 <- recode(tfr$sm4, "5 = 0; 4 = 0.2; 3 = 0.4; 2 = 0.8; 1 = 1; else = copy")
describe(tfr$SM4)


tfr$SM5 <- NA
tfr$SM5 <- recode(tfr$sm5, "5 = 0; 4 = 0.2; 3 = 0.4; 2 = 0.8; 1 = 1; else = copy")
describe(tfr$SM5)

tfr$SM6 <- NA
tfr$SM6 <- recode(tfr$sm6, "5 = 0; 4 = 0.2; 3 = 0.4; 2 = 0.8; 1 = 1; else = copy")
describe(tfr$SM6)

tfr$SM7 <- NA
tfr$SM7 <- recode(tfr$sm7, "5 = 0; 4 = 0.2; 3 = 0.4; 2 = 0.8; 1 = 1; else = copy")
describe(tfr$SM7)

tfr$SM8 <- NA
tfr$SM8 <- recode(tfr$sm8, "5 = 0; 4 = 0.2; 3 = 0.4; 2 = 0.8; 1 = 1; else = copy")
describe(tfr$SM8)


tfr$SM9 <- NA
tfr$SM9 <- recode(tfr$sm9, "5 = 0; 4 = 0.2; 3 = 0.4; 2 = 0.8; 1 = 1; else = copy")
describe(tfr$SM9)


tfr$SM10 <- NA
tfr$SM10 <- recode(tfr$sm10, "5 = 0; 4 = 0.2; 3 = 0.4; 2 = 0.8; 1 = 1; else = copy")
describe(tfr$SM10)


tfr$SM11 <- NA
tfr$SM11 <- recode(tfr$sm11, "5 = 0; 4 = 0.2; 3 = 0.4; 2 = 0.8; 1 = 1; else = copy")
describe(tfr$SM11)

tfr$SM12 <- NA
tfr$SM12 <- recode(tfr$sm12, "5 = 0; 4 = 0.2; 3 = 0.4; 2 = 0.8; 1 = 1; else = copy")
describe(tfr$SM12)



##For CM
## for CM
tfr <- rename(tfr, cm1 = CM1)
tfr <- rename(tfr, cm2 = CM2)
tfr <- rename(tfr, cm3 = CM3)
tfr <- rename(tfr, cm4 = CM4)
colnames (tfr)

tfr$CM1 <- NA
tfr$CM1 <- recode(tfr$cm1, "5 = 0; 4 = 0.2; 3 = 0.4; 2 = 0.8; 1 = 1; else = copy")
describe(tfr$CM1)


tfr$CM2 <- NA
tfr$CM2 <- recode(tfr$cm2, "5 = 0; 4 = 0.2; 3 = 0.4; 2 = 0.8; 1 = 1; else = copy")
describe(tfr$CM2)


tfr$CM3 <- NA
tfr$CM3 <- recode(tfr$cm3, "5 = 0; 4 = 0.2; 3 = 0.4; 2 = 0.8; 1 = 1; else = copy")
describe(tfr$CM3)


tfr$CM4 <- NA
tfr$CM4 <- recode(tfr$cm4, "5 = 0; 4 = 0.2; 3 = 0.4; 2 = 0.8; 1 = 1; else = copy")
describe(tfr$CM4)

######AGGREGATION#########

#build aggregated sets using the maximum rule, check for values of 0.5, missings, and for skewedness

# for W

tfr$W <- pmax(tfr$W1, tfr$W2, tfr$W3, tfr$W4, tfr$W5, na.rm=TRUE)
tfr$W

checkw <- as.numeric(tfr$W == 0.5)
sum(checkw, na.rm=TRUE)
describe(tfr$W)
sum(is.na(tfr$W))

wskew <- as.numeric(tfr$W > 0.5)
tabw <- table(wskew)
prop.table(tabw)

# for sp

tfr$W <- pmax(tfr$W1, tfr$W2, tfr$W3, tfr$W4, tfr$W5, na.rm=TRUE)
tfr$W

checkw <- as.numeric(tfr$W == 0.5)
sum(checkw, na.rm=TRUE)
describe(tfr$W)
sum(is.na(tfr$W))

wskew <- as.numeric(tfr$W > 0.5)
tabw <- table(wskew)
prop.table(tabw)

# for sp

tfr$SP <- pmax(tfr$SP1, tfr$SP2, tfr$SP3, tfr$SP4, tfr$SP5,tfr$SP6,  na.rm=TRUE)
tfr$SP

checksp <- as.numeric(tfr$SP == 0.5)
sum(checksp, na.rm=TRUE)
describe(tfr$SP)
sum(is.na(tfr$SP))

spskew <- as.numeric(tfr$SP > 0.5)
tabsp <- table(spskew)
prop.table(tabsp)

# for tp

tfr$TP <- pmax(tfr$TP1, tfr$TP2, tfr$TP3, tfr$TP4, tfr$TP5,tfr$TP6,  na.rm=TRUE)
tfr$TP

checktp <- as.numeric(tfr$TP == 0.5)
sum(checktp, na.rm=TRUE)
describe(tfr$TP)
sum(is.na(tfr$TP))

tpskew <- as.numeric(tfr$TP > 0.5)
tabtp <- table(tpskew)
prop.table(tabtp)

# for op

tfr$OP <- pmax(tfr$OP1, tfr$OP2, tfr$OP3, tfr$OP4, tfr$OP5,tfr$OP6,  na.rm=TRUE)
tfr$OP

checkop <- as.numeric(tfr$OP == 0.5)
sum(checkop, na.rm=TRUE)
describe(tfr$OP)
sum(is.na(tfr$OP))

opskew <- as.numeric(tfr$OP > 0.5)
tabop <- table(opskew)
prop.table(tabop)


# for sm

tfr$SM <- pmax(tfr$SM1, tfr$SM2, tfr$SM3, tfr$SM4, tfr$SM5,tfr$SM6, tfr$SM7, tfr$SM8, tfr$SM9, tfr$SM10, tfr$SM11,tfr$SM12, na.rm=TRUE)
tfr$SM

checksm <- as.numeric(tfr$SM == 0.5)
sum(checksm, na.rm=TRUE)
describe(tfr$SM)
sum(is.na(tfr$SM))

smskew <- as.numeric(tfr$SM > 0.5)
tabsm <- table(smskew)
prop.table(tabsm)

# for cm


tfr$CM <- pmax(tfr$CM1, tfr$CM2, tfr$CM3, tfr$CM4, na.rm = TRUE)
tfr$CM

checkcm <- as.numeric(tfr$CM == 0.5)
sum(checkcm, na.rm=TRUE)
describe(tfr$CM)
sum(is.na(tfr$CM))

cmskew <- as.numeric(tfr$CM > 0.5)
tabcm <- table(cmskew)
prop.table(tabcm)


#exclude cases with set membership 0.5

tfr = tfr[tfr$W != 0.5,]
tfr = tfr[tfr$OP != 0.5,] 
tfr = tfr[tfr$SP != 0.5,]
tfr = tfr[tfr$TP != 0.5,]
tfr = tfr[tfr$SM != 0.5,]
tfr = tfr[tfr$CM != 0.5,]
tfr

describe(tfr)


tff_d1_c2 <- subset(tfr, select = c("W", "SP", "TP", "OP", "SM", "CM"))
write.csv2(tff_d1_c2, file = "tff_d1_c2.csv")


###################ANALYSIS OF NECESSITY############

tff <- read.csv("tff_d1_c2.csv", row.names=1, header = TRUE, sep = ";",  dec = ",")
describe(tff)

#analysis of necessity for W and w  

superSubset(tff, outcome = "W", 
            conditions = "OP, SP, TP, SM, CM", 
            incl.cut = 0.9, cov.cut = 0.6)

superSubset(tff, outcome = "~W", 
            conditions = "OP, SP, TP, SM, CM", 
            incl.cut = 0.9, cov.cut = 0.6)




#calculate Z-Test for necessary conditions. I apply a RoN threshold of 0.55, so only for NCs that meet this threshold

nw <- 1-tff$W

# for NC9 OP + TP

NC9W <- compute("OP+TP", data=tff)
obsNC9W <- as.numeric(NC9W >= tff$W)

zNC9W <- round(((((sum(obsNC9W == 1))/1004) - 0.8)-(1/(2*1004))) / (sqrt((0.8*(1-0.8))/1004)), digits = 1)
zNC9W

pzNC9W <- 2*pnorm(-abs(zNC9W))
pzNC9W

# for NC11 OP+SP

NC11W <- compute("OP+SP", data=tff)
obsNC11W <- as.numeric(NC11W >= tff$W)

zNC11W <- round(((((sum(obsNC11W == 1))/1004) - 0.8)-(1/(2*1004))) / (sqrt((0.8*(1-0.8))/1004)), digits = 1)
zNC11W

pzNC11W <- 2*pnorm(-abs(zNC11W))
pzNC11W

# for NC16 SP+TP+SM

NC16W <- compute("SP+TP+SM", data=tff)
obsNC16W <- as.numeric(NC16W >= tff$W)

zNC16W <-round( ((((sum(obsNC16W == 1))/1004) - 0.8)-(1/(2*1004))) / (sqrt((0.8*(1-0.8))/1004)), digits=1)
zNC16W

pzNC16W <- 2*pnorm(-abs(zNC16W))
pzNC16W

# for NC19  op+SP+SM

NC19W <- compute("op+SP+SM", data=tff)
obsNC19W <- as.numeric(NC19W >= tff$W)

zNC19W <-round( ((((sum(obsNC19W == 1))/1004) - 0.8)-(1/(2*1004))) / (sqrt((0.8*(1-0.8))/1004)), digits=1)
zNC19W

pzNC19W <- 2*pnorm(-abs(zNC19W))
pzNC19W

# for NC22  OP+SM+CM

NC22W <- compute("OP+SM+CM", data=tff)
obsNC22W <- as.numeric(NC22W >= tff$W)

zNC22W <-round( ((((sum(obsNC22W == 1))/1004) - 0.8)-(1/(2*1004))) / (sqrt((0.8*(1-0.8))/1004)), digits=1)
zNC22W

pzNC22W <- 2*pnorm(-abs(zNC22W))
pzNC22W




#Test for skewness of necessary conditions. It is tested whether the proportion of cases in all but one 
#specific area of the XY plot is lower than ten per cent.
##First, visual check with XY Plots; then calculate number of cases outside each area; 
#and caluclate percentage outside area; 
#there are a lot of cases on the diagonal, which are attributed to two area. This is why 
#the percentages don't add up to 100.


##for TP+OP

A1NC9W <- round(sum(as.numeric((NC9W <= tff$W) & (NC9W >= nw)))/(1004/100), digits = 1)

A2NC9W <- round(sum(as.numeric((NC9W >= tff$W)&(NC9W >= nw)))/(1004/100), digits = 1)

A3NC9W <- round(sum(as.numeric((NC9W >= tff$W)&(NC9W <= nw)))/(1004/100), digits=1)

A4NC9W <- round(sum(as.numeric((NC9W <= tff$W)&(NC9W <= nw)))/(1004/100), digits = 1)

sum(A1NC9W, A2NC9W, A3NC9W, A4NC9W)

par(mfrow=c(3, 2))
plot(tff$W, NC9W, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main='Skewness test for TP + OP',
     xlab='TP + OP',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1NC9W)
text(x=0.8, y=0.5, labels = A2NC9W)
text(x=0.5, y=0.2, labels = A3NC9W)
text(x=0.2, y=0.5, labels = A4NC9W)

# for NC11 SP + OP

A1NC11W <- round(sum(as.numeric((NC11W <= tff$W) & (NC11W >= nw)))/(1004/100), digits = 1)

A2NC11W <- round(sum(as.numeric((NC11W >= tff$W)&(NC11W >= nw)))/(1004/100), digits = 1)

A3NC11W <- round(sum(as.numeric((NC11W >= tff$W)&(NC11W <= nw)))/(1004/100), digits=1)

A4NC11W <- round(sum(as.numeric((NC11W <= tff$W)&(NC11W <= nw)))/(1004/100), digits = 1)

sum(A1NC11W, A2NC11W, A3NC11W, A4NC11W)

plot(tff$W, NC11W, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main='Skewness test for SP + OP',
     xlab='SP + OP',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1NC11W)
text(x=0.8, y=0.5, labels = A2NC11W)
text(x=0.5, y=0.2, labels = A3NC11W)
text(x=0.2, y=0.5, labels = A4NC11W)

##for NC 16 SP+TP+SM

A1NC16W <- round(sum(as.numeric((NC16W <= tff$W) & (NC16W >= nw)))/(1004/100), digits = 1)

A2NC16W <- round(sum(as.numeric((NC16W >= tff$W)&(NC16W >= nw)))/(1004/100), digits = 1)

A3NC16W <- round(sum(as.numeric((NC16W >= tff$W)&(NC16W <= nw)))/(1004/100), digits=1)

A4NC16W <- round(sum(as.numeric((NC16W <= tff$W)&(NC16W <= nw)))/(1004/100), digits = 1)

sum(A1NC16W, A2NC16W, A3NC16W, A4NC16W)

plot(tff$W, NC16W, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main='Skewness test for SP+TP+SM',
     xlab=' SP+TP+SM ',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1NC16W)
text(x=0.8, y=0.5, labels = A2NC16W)
text(x=0.5, y=0.2, labels = A3NC16W)
text(x=0.2, y=0.5, labels = A4NC16W)


# for NC19  op+SP+SM

A1NC19W <- round(sum(as.numeric((NC19W <= tff$W) & (NC19W >= nw)))/(1004/100), digits = 1)

A2NC19W <- round(sum(as.numeric((NC19W >= tff$W)&(NC19W >= nw)))/(1004/100), digits = 1)

A3NC19W <- round(sum(as.numeric((NC19W >= tff$W)&(NC19W <= nw)))/(1004/100), digits=1)

A4NC19W <- round(sum(as.numeric((NC19W <= tff$W)&(NC19W <= nw)))/(1004/100), digits = 1)

sum(A1NC19W, A2NC19W, A3NC19W, A4NC19W)

plot(tff$W, NC19W, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main='Skewness test for op+SP+SM',
     xlab='op+SP+SM',
     ylab='W')
     abline(0,1, col="black")
     abline(1, -1, col="black")
     text(x=0.5, y=0.8, labels = A1NC19W)
     text(x=0.8, y=0.5, labels = A2NC19W)
     text(x=0.5, y=0.2, labels = A3NC19W)
     text(x=0.2, y=0.5, labels = A4NC19W)
     

     # for NC22  OP+SM+CM
     
     A1NC22W <- round(sum(as.numeric((NC22W <= tff$W) & (NC22W >= nw)))/(1004/100), digits = 1)
     
     A2NC22W <- round(sum(as.numeric((NC22W >= tff$W)&(NC22W >= nw)))/(1004/100), digits = 1)
     
     A3NC22W <- round(sum(as.numeric((NC22W >= tff$W)&(NC22W <= nw)))/(1004/100), digits=1)
     
     A4NC22W <- round(sum(as.numeric((NC22W <= tff$W)&(NC22W <= nw)))/(1004/100), digits = 1)
     
     sum(A1NC22W, A2NC22W, A3NC22W, A4NC22W)
     
     plot(tff$W, NC22W, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
          xaxs="i", yaxs="i",
          main='Skewness test for OP+SM+CM ',
          xlab=' OP+SM+CM ',
          ylab='W')
     abline(0,1, col="black")
     abline(1, -1, col="black")
     text(x=0.5, y=0.8, labels = A1NC22W)
     text(x=0.8, y=0.5, labels = A2NC22W)
     text(x=0.5, y=0.2, labels = A3NC22W)
     text(x=0.2, y=0.5, labels = A4NC22W)

#####check accuracy of untenable assumptions####
     
#as the untenable assumptions for the analysis of sufficiency for W were set sp*op + tp*op, it was assumed that the two necessary conditions are connectd with the logical AND. However, with imperfect subset relations this is not always the case. 
#Check if it is the case:
     
     necall <- compute("(TP+OP)*(SP+OP)", data=tff)
     
     pof(necall, W, tff, relation = "nec")
     
# the condition has consistency necessity of .887. Hence, the untenable assumptions are slightly overly restrictive.
     
           
#################ANALYSIS OF SUFFICIENCY FOR W##################
# 


tff <- read.csv("tff_d1_c2.csv", row.names=1, header = TRUE, sep = ";",  dec = ",")
describe(tff)

nw <- 1-tff$W

#analysis of sufficiency for W
# calculate untenable assumptions

nec <- deMorgan("(TP + OP)*(SP + OP)", snames = "TP, OP, SP")
nec

#### W1: frequency threshold 5, raw consistency threshold 0.9468 ####

ttW1 <- truthTable(tff, outcome="W", 
                        conditions = "SP, TP, OP, SM, CM", 
                        incl.cut=0.9468, n.cut=5, sort.by="incl, n", decreasing=TRUE, 
                        complete=FALSE, show.cases=FALSE)
ttW1


colnames(tff)

#parsimonious solution
psW1 <- eqmcc(ttW1, include="?", details=TRUE, show.cases=FALSE,  row.dom=TRUE, all.sol=FALSE)

psW1

#exclude contradictory assumptions & rows: tp*op + sp*op (contradict statement of necessity)

ttW1new <- ttW1
nec
rows <- findRows("tp*op + op*sp", ttW1new, remainders = FALSE)
rows

ttW1new$tt[as.character(rows), "OUT"] <- 0
ttW1new


#conservative solution
csW1 <- eqmcc(ttW1new, details=TRUE, PRI=TRUE, show.cases=FALSE, rowdom=TRUE, use.tilde=FALSE)
csW1

#enhanced parsimonious solution (all simplifying assumptions)

epsW1 <- eqmcc(ttW1new, include="?", details=TRUE, PRI=TRUE, show.cases=FALSE, 
              rowdom=TRUE, use.tilde=FALSE)
epsW1


#Z-Test for enhanced parsimonious solution, 0.8 (almost always sufficient)
# Calculate membership in solution term  sp*OP*CM + tp*OP*CM + SP*TP*op*CM + (SP*TP*SM*CM)
epsolW1 <- compute("sp*OP*CM + tp*OP*CM + SP*TP*op*CM + SP*TP*SM*CM", data=tff)

obsepsolW1 <- as.numeric(epsolW1 <= tff$W)
zepsolW1 <- round(((((sum(obsepsolW1 == 1))/1004) - 0.8)-(1/(2*1004))) / (sqrt((0.8*(1-0.8))/1004)), digits=1)
zepsolW1
pzepsolW1 <- 2*pnorm(-abs(zepsolW1))
pzepsolW1


#intermediate solution with directional expectations SM -> W, CM -> W
isW1 <- eqmcc(ttW1new, details = TRUE, include = "?",  use.tilde = FALSE, PRI = TRUE,
              row.dom = TRUE, show.cases = FALSE, 
              dir.exp = "-, -, -, 1, 1")

isW1



#Test for skewness of sufficient conditions (enhanced parsimonious solution)
# calculate membership in paths of M2: sp*OP*CM + tp*OP*CM + SP*TP*op*CM + (SP*TP*SM*CM)
path1W1 <- compute("sp*OP*CM", data=tff)
path2W1 <- compute("tp*OP*CM", data=tff)
path3W1 <- compute("SP*TP*op*CM", data=tff)
path4W1 <- compute("SP*TP*SM*CM", data=tff)

A1path1W1 <- round(sum(as.numeric((path1W1 <= tff$W) & (path1W1 >= nw)))/(1004/100), digits = 1)

A2path1W1 <- round(sum(as.numeric((path1W1 >= tff$W)&(path1W1 >= nw)))/(1004/100), digits = 1)

A3path1W1 <- round(sum(as.numeric((path1W1 >= tff$W)&(path1W1 <= nw)))/(1004/100), digits=1)

A4path1W1 <- round(sum(as.numeric((path1W1 <= tff$W)&(path1W1 <= nw)))/(1004/100), digits = 1)


par(mfrow=c(2, 3))
plot(tff$W, path1W1, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 1',
     xlab= 'Path 1',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path1W1)
text(x=0.8, y=0.5, labels = A2path1W1)
text(x=0.5, y=0.2, labels = A3path1W1)
text(x=0.2, y=0.5, labels = A4path1W1)



A1path2W1 <- round(sum(as.numeric((path2W1 <= tff$W) & (path2W1 >= nw)))/(1004/100), digits = 1)

A2path2W1 <- round(sum(as.numeric((path2W1 >= tff$W)&(path2W1 >= nw)))/(1004/100), digits = 1)

A3path2W1 <- round(sum(as.numeric((path2W1 >= tff$W)&(path2W1 <= nw)))/(1004/100), digits=1)

A4path2W1 <- round(sum(as.numeric((path2W1 <= tff$W)&(path2W1 <= nw)))/(1004/100), digits = 1)

plot(tff$W, path2W1, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 2',
     xlab='Path 2',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path2W1)
text(x=0.8, y=0.5, labels = A2path2W1)
text(x=0.5, y=0.2, labels = A3path2W1)
text(x=0.2, y=0.5, labels = A4path2W1)

A1path3W1 <- round(sum(as.numeric((path3W1 <= tff$W) & (path3W1 >= nw)))/(1004/100), digits = 1)

A2path3W1 <- round(sum(as.numeric((path3W1 >= tff$W)&(path3W1 >= nw)))/(1004/100), digits = 1)

A3path3W1 <- round(sum(as.numeric((path3W1 >= tff$W)&(path3W1 <= nw)))/(1004/100), digits=1)

A4path3W1 <- round(sum(as.numeric((path3W1 <= tff$W)&(path3W1 <= nw)))/(1004/100), digits = 1)

plot(tff$W, path3W1, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main='Skewness path 3',
     xlab='Path 3',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path3W1)
text(x=0.8, y=0.5, labels = A2path3W1)
text(x=0.5, y=0.2, labels = A3path3W1)
text(x=0.2, y=0.5, labels = A4path3W1)

A1path4W1 <- round(sum(as.numeric((path4W1 <= tff$W) & (path4W1 >= nw)))/(1004/100), digits = 1)

A2path4W1 <- round(sum(as.numeric((path4W1 >= tff$W)&(path4W1 >= nw)))/(1004/100), digits = 1)

A3path4W1 <- round(sum(as.numeric((path4W1 >= tff$W)&(path4W1 <= nw)))/(1004/100), digits=1)

A4path4W1 <- round(sum(as.numeric((path4W1 <= tff$W)&(path4W1 <= nw)))/(1004/100), digits = 1)

plot(tff$W, path4W1, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main='Skewness path 4',
     xlab='Path 4',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path4W1)
text(x=0.8, y=0.5, labels = A2path4W1)
text(x=0.5, y=0.2, labels = A3path4W1)
text(x=0.2, y=0.5, labels = A4path4W1)

A1epsolW1 <- round(sum(as.numeric((epsolW1 <= tff$W) & (epsolW1 >= nw)))/(1004/100), digits = 1)

A2epsolW1 <- round(sum(as.numeric((epsolW1 >= tff$W)&(epsolW1 >= nw)))/(1004/100), digits = 1)

A3epsolW1 <- round(sum(as.numeric((epsolW1 >= tff$W)&(epsolW1 <= nw)))/(1004/100), digits=1)

A4epsolW1 <- round(sum(as.numeric((epsolW1 <= tff$W)&(epsolW1 <= nw)))/(1004/100), digits = 1)

isW1
plot(tff$W, epsolW1, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main=' Skewness eps',
     xlab='Solution term',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1epsolW1)
text(x=0.8, y=0.5, labels = A2epsolW1)
text(x=0.5, y=0.2, labels = A3epsolW1)
text(x=0.2, y=0.5, labels = A4epsolW1)

#####W2  frequency threshold 5, raw consistency threshold 0.9405 #####


ttW2 <- truthTable(tff, outcome="W", 
                   conditions = "SP, TP, OP, SM, CM", 
                   incl.cut=0.9405, n.cut=5, sort.by="incl, n", decreasing=TRUE, 
                   complete=FALSE, show.cases=FALSE)
ttW2


colnames(tff)

#parsimonious solution
psW2 <- eqmcc(ttW2, include="?", details=TRUE, show.cases=FALSE,  row.dom=TRUE, all.sol=FALSE)

psW2

#exclude contradictory assumptions & rows: tp*op + sp*op (contradict statement of necessity)

ttW2new <- ttW2
nec
rows <- findRows("tp*op + op*sp", ttW2new, remainders = FALSE)
rows

ttW2new$tt[as.character(rows), "OUT"] <- 0
ttW2new


#conservative solution
csW2 <- eqmcc(ttW2new, details=TRUE, PRI=TRUE, show.cases=FALSE, rowdom=TRUE, use.tilde=FALSE)
csW2

#enhanced parsimonious solution (all simplifying assumptions)

epsW2 <- eqmcc(ttW2new, include="?", details=TRUE, PRI=TRUE, show.cases=FALSE, 
               rowdom=TRUE, use.tilde=FALSE)
epsW2


#Z-Test for enhanced parsimonious solution, 0.8 (almost always sufficient)
# Calculate membership in solution term  sp*OP*CM + tp*OP*CM + SP*TP*op*CM + sp*TP*OP*SM + SP*TP*op*SM + SP*TP*SM*CM
epsolW2 <- compute("sp*OP*CM + tp*OP*CM + SP*TP*op*CM + sp*TP*OP*SM + SP*TP*op*SM + SP*TP*SM*CM ", data=tff)

obsepsolW2 <- as.numeric(epsolW2 <= tff$W)
zepsolW2 <- round(((((sum(obsepsolW2 == 1))/1004) - 0.8)-(1/(2*1004))) / (sqrt((0.8*(1-0.8))/1004)), digits=1)
zepsolW2
pzepsolW2 <- 2*pnorm(-abs(zepsolW2))
pzepsolW2


#intermediate solution with directional expectations SM -> W, CM -> W
isW2 <- eqmcc(ttW2new, details = TRUE, include = "?",  use.tilde = FALSE, PRI = TRUE,
              row.dom = TRUE, show.cases = FALSE, 
              dir.exp = "-, -, -, 1, 1")

isW2




#Test for skewness of sufficient conditions (enhanced parsimonious solution)
# calculate membership in paths of M2: sp*OP*CM + tp*OP*CM + SP*TP*op*CM + sp*TP*OP*SM + SP*TP*op*SM + SP*TP*SM*CM
path1W2 <- compute("sp*OP*CM ", data=tff)
path2W2 <- compute("tp*OP*CM ", data=tff)
path3W2 <- compute("SP*TP*op*CM ", data=tff)
path4W2 <- compute("sp*TP*OP*SM ", data=tff)
path5W2 <- compute("SP*TP*op*SM ", data=tff)
path6W2 <- compute("SP*TP*SM*CM ", data=tff)

A1path1W2 <- round(sum(as.numeric((path1W2 <= tff$W) & (path1W2 >= nw)))/(1004/100), digits = 1)

A2path1W2 <- round(sum(as.numeric((path1W2 >= tff$W)&(path1W2 >= nw)))/(1004/100), digits = 1)

A3path1W2 <- round(sum(as.numeric((path1W2 >= tff$W)&(path1W2 <= nw)))/(1004/100), digits=1)

A4path1W2 <- round(sum(as.numeric((path1W2 <= tff$W)&(path1W2 <= nw)))/(1004/100), digits = 1)


par(mfrow=c(2, 3))
plot(tff$W, path1W2, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 1',
     xlab= 'Path 1',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path1W2)
text(x=0.8, y=0.5, labels = A2path1W2)
text(x=0.5, y=0.2, labels = A3path1W2)
text(x=0.2, y=0.5, labels = A4path1W2)



A1path2W2 <- round(sum(as.numeric((path2W2 <= tff$W) & (path2W2 >= nw)))/(1004/100), digits = 1)

A2path2W2 <- round(sum(as.numeric((path2W2 >= tff$W)&(path2W2 >= nw)))/(1004/100), digits = 1)

A3path2W2 <- round(sum(as.numeric((path2W2 >= tff$W)&(path2W2 <= nw)))/(1004/100), digits=1)

A4path2W2 <- round(sum(as.numeric((path2W2 <= tff$W)&(path2W2 <= nw)))/(1004/100), digits = 1)

plot(tff$W, path2W2, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 2',
     xlab='Path 2',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path2W2)
text(x=0.8, y=0.5, labels = A2path2W2)
text(x=0.5, y=0.2, labels = A3path2W2)
text(x=0.2, y=0.5, labels = A4path2W2)

A1path3W2 <- round(sum(as.numeric((path3W2 <= tff$W) & (path3W2 >= nw)))/(1004/100), digits = 1)

A2path3W2 <- round(sum(as.numeric((path3W2 >= tff$W)&(path3W2 >= nw)))/(1004/100), digits = 1)

A3path3W2 <- round(sum(as.numeric((path3W2 >= tff$W)&(path3W2 <= nw)))/(1004/100), digits=1)

A4path3W2 <- round(sum(as.numeric((path3W2 <= tff$W)&(path3W2 <= nw)))/(1004/100), digits = 1)

plot(tff$W, path3W2, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main='Skewness path 3',
     xlab='Path 3',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path3W2)
text(x=0.8, y=0.5, labels = A2path3W2)
text(x=0.5, y=0.2, labels = A3path3W2)
text(x=0.2, y=0.5, labels = A4path3W2)

A1path4W2 <- round(sum(as.numeric((path4W2 <= tff$W) & (path4W2 >= nw)))/(1004/100), digits = 1)

A2path4W2 <- round(sum(as.numeric((path4W2 >= tff$W)&(path4W2 >= nw)))/(1004/100), digits = 1)

A3path4W2 <- round(sum(as.numeric((path4W2 >= tff$W)&(path4W2 <= nw)))/(1004/100), digits=1)

A4path4W2 <- round(sum(as.numeric((path4W2 <= tff$W)&(path4W2 <= nw)))/(1004/100), digits = 1)

plot(tff$W, path4W2, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main='Skewness path 4',
     xlab='Path 4',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path4W2)
text(x=0.8, y=0.5, labels = A2path4W2)
text(x=0.5, y=0.2, labels = A3path4W2)
text(x=0.2, y=0.5, labels = A4path4W2)

A1path5W2 <- round(sum(as.numeric((path5W2 <= tff$W) & (path5W2 >= nw)))/(1004/100), digits = 1)

A2path5W2 <- round(sum(as.numeric((path5W2 >= tff$W)&(path5W2 >= nw)))/(1004/100), digits = 1)

A3path5W2 <- round(sum(as.numeric((path5W2 >= tff$W)&(path5W2 <= nw)))/(1004/100), digits=1)

A4path5W2 <- round(sum(as.numeric((path5W2 <= tff$W)&(path5W2 <= nw)))/(1004/100), digits = 1)


par(mfrow=c(2, 3))
plot(tff$W, path5W2, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 5',
     xlab= 'Path 5',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path5W2)
text(x=0.8, y=0.5, labels = A2path5W2)
text(x=0.5, y=0.2, labels = A3path5W2)
text(x=0.2, y=0.5, labels = A4path5W2)



A1path6W2 <- round(sum(as.numeric((path6W2 <= tff$W) & (path6W2 >= nw)))/(1004/100), digits = 1)

A2path6W2 <- round(sum(as.numeric((path6W2 >= tff$W)&(path6W2 >= nw)))/(1004/100), digits = 1)

A3path6W2 <- round(sum(as.numeric((path6W2 >= tff$W)&(path6W2 <= nw)))/(1004/100), digits=1)

A4path6W2 <- round(sum(as.numeric((path6W2 <= tff$W)&(path6W2 <= nw)))/(1004/100), digits = 1)

plot(tff$W, path6W2, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 6',
     xlab='Path 6',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path6W2)
text(x=0.8, y=0.5, labels = A2path6W2)
text(x=0.5, y=0.2, labels = A3path6W2)
text(x=0.2, y=0.5, labels = A4path6W2)



A1epsolW2 <- round(sum(as.numeric((epsolW2 <= tff$W) & (epsolW2 >= nw)))/(1004/100), digits = 1)

A2epsolW2 <- round(sum(as.numeric((epsolW2 >= tff$W)&(epsolW2 >= nw)))/(1004/100), digits = 1)

A3epsolW2 <- round(sum(as.numeric((epsolW2 >= tff$W)&(epsolW2 <= nw)))/(1004/100), digits=1)

A4epsolW2 <- round(sum(as.numeric((epsolW2 <= tff$W)&(epsolW2 <= nw)))/(1004/100), digits = 1)

isW2
plot(tff$W, epsolW2, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main=' Skewness eps',
     xlab='Solution term',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1epsolW2)
text(x=0.8, y=0.5, labels = A2epsolW2)
text(x=0.5, y=0.2, labels = A3epsolW2)
text(x=0.2, y=0.5, labels = A4epsolW2)


####W3: frequency threshold 5, raw consistency threshold 0.9395#######

ttW3 <- truthTable(tff, outcome="W", 
                   conditions = "SP, TP, OP, SM, CM", 
                   incl.cut=0.9395, n.cut=5, sort.by="incl, n", decreasing=TRUE, 
                   complete=FALSE, show.cases=FALSE)
ttW3


#parsimonious solution
psW3 <- eqmcc(ttW3, include="?", details=TRUE, show.cases=FALSE,  row.dom=TRUE, all.sol=FALSE)

psW3

#exclude contradictory assumptions & rows: tp*op + sp*op (contradict statement of necessity)

ttW3new <- ttW3
nec
rows <- findRows("tp*op + op*sp", ttW3new, remainders = FALSE)
rows

ttW3new$tt[as.character(rows), "OUT"] <- 0
ttW3new


#conservative solution
csW3 <- eqmcc(ttW3new, details=TRUE, PRI=TRUE, show.cases=FALSE, rowdom=TRUE, use.tilde=FALSE)
csW3

#enhanced parsimonious solution (all simplifying assumptions)

epsW3 <- eqmcc(ttW3new, include="?", details=TRUE, PRI=TRUE, show.cases=FALSE, 
               rowdom=TRUE, use.tilde=FALSE)
epsW3

#Z-Test for enhanced parsimonious solution, 0.8 (almost always sufficient)
# Calculate membership in solution term OP*CM + SP*TP*CM + sp*TP*OP*SM + SP*TP*op*SM
epsolW3 <- compute("OP*CM + SP*TP*CM + sp*TP*OP*SM + SP*TP*op*SM", data=tff)

obsepsolW3 <- as.numeric(epsolW3 <= tff$W)
zepsolW3 <- round(((((sum(obsepsolW3 == 1))/1004) - 0.8)-(1/(2*1004))) / (sqrt((0.8*(1-0.8))/1004)), digits=1)
zepsolW3
pzepsolW3 <- 2*pnorm(-abs(zepsolW3))
pzepsolW3


#intermediate solution with directional expectations SM -> W, CM -> W
isW3 <- eqmcc(ttW3new, details = TRUE, include = "?",  use.tilde = FALSE, PRI = TRUE,
              row.dom = TRUE, show.cases = FALSE, 
              dir.exp = "-, -, -, 1, 1")

isW3

# identify simplifying assumptions and easy counterfactuals

SAW3 <- psW3$SA
SAW3

ECW3 <- isW3$i.sol$C1P1$EC
ECW3


#Test for skewness of sufficient conditions (enhanced parsimonious solution)
# calculate membership in paths of M1: OP*CM + SP*TP*CM + sp*TP*OP*SM + SP*TP*op*SM

path1W3 <- compute("OP*CM", data=tff)
path2W3 <- compute("SP*TP*CM ", data=tff)
path3W3 <- compute("sp*TP*OP*SM ", data=tff)
path4W3 <- compute("SP*TP*op*SM ", data=tff)

A1path1W3 <- round(sum(as.numeric((path1W3 <= tff$W) & (path1W3 >= nw)))/(1004/100), digits = 1)

A2path1W3 <- round(sum(as.numeric((path1W3 >= tff$W)&(path1W3 >= nw)))/(1004/100), digits = 1)

A3path1W3 <- round(sum(as.numeric((path1W3 >= tff$W)&(path1W3 <= nw)))/(1004/100), digits=1)

A4path1W3 <- round(sum(as.numeric((path1W3 <= tff$W)&(path1W3 <= nw)))/(1004/100), digits = 1)


par(mfrow=c(2, 3))
plot(tff$W, path1W3, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 1',
     xlab= 'Path 1',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path1W3)
text(x=0.8, y=0.5, labels = A2path1W3)
text(x=0.5, y=0.2, labels = A3path1W3)
text(x=0.2, y=0.5, labels = A4path1W3)



A1path2W3 <- round(sum(as.numeric((path2W3 <= tff$W) & (path2W3 >= nw)))/(1004/100), digits = 1)

A2path2W3 <- round(sum(as.numeric((path2W3 >= tff$W)&(path2W3 >= nw)))/(1004/100), digits = 1)

A3path2W3 <- round(sum(as.numeric((path2W3 >= tff$W)&(path2W3 <= nw)))/(1004/100), digits=1)

A4path2W3 <- round(sum(as.numeric((path2W3 <= tff$W)&(path2W3 <= nw)))/(1004/100), digits = 1)

plot(tff$W, path2W3, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 2',
     xlab='Path 2',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path2W3)
text(x=0.8, y=0.5, labels = A2path2W3)
text(x=0.5, y=0.2, labels = A3path2W3)
text(x=0.2, y=0.5, labels = A4path2W3)

A1path3W3 <- round(sum(as.numeric((path3W3 <= tff$W) & (path3W3 >= nw)))/(1004/100), digits = 1)

A2path3W3 <- round(sum(as.numeric((path3W3 >= tff$W)&(path3W3 >= nw)))/(1004/100), digits = 1)

A3path3W3 <- round(sum(as.numeric((path3W3 >= tff$W)&(path3W3 <= nw)))/(1004/100), digits=1)

A4path3W3 <- round(sum(as.numeric((path3W3 <= tff$W)&(path3W3 <= nw)))/(1004/100), digits = 1)

plot(tff$W, path3W3, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main='Skewness path 3',
     xlab='Path 3',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path3W3)
text(x=0.8, y=0.5, labels = A2path3W3)
text(x=0.5, y=0.2, labels = A3path3W3)
text(x=0.2, y=0.5, labels = A4path3W3)

A1path4W3 <- round(sum(as.numeric((path4W3 <= tff$W) & (path4W3 >= nw)))/(1004/100), digits = 1)

A2path4W3 <- round(sum(as.numeric((path4W3 >= tff$W)&(path4W3 >= nw)))/(1004/100), digits = 1)

A3path4W3 <- round(sum(as.numeric((path4W3 >= tff$W)&(path4W3 <= nw)))/(1004/100), digits=1)

A4path4W3 <- round(sum(as.numeric((path4W3 <= tff$W)&(path4W3 <= nw)))/(1004/100), digits = 1)

plot(tff$W, path4W3, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main='Skewness path 4',
     xlab='Path 4',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path4W3)
text(x=0.8, y=0.5, labels = A2path4W3)
text(x=0.5, y=0.2, labels = A3path4W3)
text(x=0.2, y=0.5, labels = A4path4W3)

A1epsolW3 <- round(sum(as.numeric((epsolW3 <= tff$W) & (epsolW3 >= nw)))/(1004/100), digits = 1)

A2epsolW3 <- round(sum(as.numeric((epsolW3 >= tff$W)&(epsolW3 >= nw)))/(1004/100), digits = 1)

A3epsolW3 <- round(sum(as.numeric((epsolW3 >= tff$W)&(epsolW3 <= nw)))/(1004/100), digits=1)

A4epsolW3 <- round(sum(as.numeric((epsolW3 <= tff$W)&(epsolW3 <= nw)))/(1004/100), digits = 1)

isW3
plot(tff$W, epsolW3, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main=' Skewness eps',
     xlab='Solution term',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1epsolW3)
text(x=0.8, y=0.5, labels = A2epsolW3)
text(x=0.5, y=0.2, labels = A3epsolW3)
text(x=0.2, y=0.5, labels = A4epsolW3)

#Z-Test for intermediate solution, 0.8 (almost always sufficient)
# Calculate membership in solution term  OP*CM + SP*TP*CM + sp*TP*OP*SM + SP*TP*op*SM
isolW3 <- compute("OP*CM + SP*TP*CM + sp*TP*OP*SM + SP*TP*op*SM", data=tff)

obsisolW3 <- as.numeric(isolW3 <= tff$W)
zisolW3 <- round(((((sum(obsisolW3 == 1))/1004) - 0.8)-(1/(2*1004))) / (sqrt((0.8*(1-0.8))/1004)), digits=1)
zisolW3
pzisolW3 <- 2*pnorm(-abs(zisolW3))
pzisolW3

# for path 1
# Calculate membership in path  OP*CM
path1isolW3 <- compute("OP*CM ", data=tff)

obspath1isolW3 <- as.numeric(path1isolW3 <= tff$W)
zpath1isolW3 <- round(((((sum(obspath1isolW3 == 1))/1004) - 0.8)-(1/(2*1004))) / (sqrt((0.8*(1-0.8))/1004)), digits=1)
zpath1isolW3
pzpath1isolW3 <- 2*pnorm(-abs(zpath1isolW3))
pzpath1isolW3

#for path 2
# Calculate membership in path  SP*TP*CM
path2isolW3 <- compute("SP*TP*CM ", data=tff)

obspath2isolW3 <- as.numeric(path2isolW3 <= tff$W)
zpath2isolW3 <- round(((((sum(obspath2isolW3 == 1))/1004) - 0.8)-(1/(2*1004))) / (sqrt((0.8*(1-0.8))/1004)), digits=1)
zpath2isolW3
pzpath2isolW3 <- 2*pnorm(-abs(zpath2isolW3))
pzpath2isolW3

# for path 3
# Calculate membership in path  sp*TP*OP*SM
path3isolW3 <- compute("sp*TP*OP*SM ", data=tff)

obspath3isolW3 <- as.numeric(path3isolW3 <= tff$W)
zpath3isolW3 <- round(((((sum(obspath3isolW3 == 1))/1004) - 0.8)-(1/(2*1004))) / (sqrt((0.8*(1-0.8))/1004)), digits=1)
zpath3isolW3
pzpath3isolW3 <- 2*pnorm(-abs(zpath3isolW3))
pzpath3isolW3

# for path 4
# Calculate membership in path  SP*TP*op*SM
path4isolW3 <- compute("SP*TP*op*SM ", data=tff)

obspath4isolW3 <- as.numeric(path4isolW3 <= tff$W)
zpath4isolW3 <- round(((((sum(obspath4isolW3 == 1))/1004) - 0.8)-(1/(2*1004))) / (sqrt((0.8*(1-0.8))/1004)), digits=1)
zpath4isolW3
pzpath4isolW3 <- 2*pnorm(-abs(zpath4isolW3))
pzpath4isolW3


####W4 frequency threshold 10, raw consistency threshold 0.9405####


ttW4 <- truthTable(tff, outcome="W", 
                   conditions = "SP, TP, OP, SM, CM", 
                   incl.cut=0.9405, n.cut=10, sort.by="incl, n", decreasing=TRUE, 
                   complete=FALSE, show.cases=FALSE)
ttW4


colnames(tff)

#parsimonious solution
psW4 <- eqmcc(ttW4, include="?", details=TRUE, show.cases=FALSE,  row.dom=TRUE, all.sol=FALSE)

psW4

#exclude contradictory assumptions & rows: tp*op + sp*op (contradict statement of necessity)

ttW4new <- ttW4
nec
rows <- findRows("tp*op + op*sp", ttW4new, remainders = FALSE)
rows

ttW4new$tt[as.character(rows), "OUT"] <- 0
ttW4new


#conservative solution
csW4 <- eqmcc(ttW4new, details=TRUE, PRI=TRUE, show.cases=FALSE, rowdom=TRUE, use.tilde=FALSE)
csW4

#enhanced parsimonious solution (all simplifying assumptions)

epsW4 <- eqmcc(ttW4new, include="?", details=TRUE, PRI=TRUE, show.cases=FALSE, 
               rowdom=TRUE, use.tilde=FALSE)
epsW4


#Z-Test for enhanced parsimonious solution, 0.8 (almost always sufficient)
# Calculate membership in solution term OP*SM*CM + sp*TP*OP*SM + tp*OP*CM
epsolW4 <- compute("OP*SM*CM + sp*TP*OP*SM + tp*OP*CM", data=tff)

obsepsolW4 <- as.numeric(epsolW4 <= tff$W)
zepsolW4 <- round(((((sum(obsepsolW4 == 1))/1004) - 0.8)-(1/(2*1004))) / (sqrt((0.8*(1-0.8))/1004)), digits=1)
zepsolW4
pzepsolW4 <- 2*pnorm(-abs(zepsolW4))
pzepsolW4


#intermediate solution with directional expectations SM -> W, CM -> W
isW4 <- eqmcc(ttW4new, details = TRUE, include = "?",  use.tilde = FALSE, PRI = TRUE,
              row.dom = TRUE, show.cases = FALSE, 
              dir.exp = "-, -, -, 1, 1")

isW4

#Test for skewness of sufficient conditions (enhanced parsimonious solution)
# calculate membership in paths of M2: OP*SM*CM + sp*TP*OP*SM + tp*OP*CM

path1W4 <- compute("OP*SM*CM ", data=tff)
path2W4 <- compute("sp*TP*OP*SM ", data=tff)
path3W4 <- compute("tp*OP*CM ", data=tff)

A1path1W4 <- round(sum(as.numeric((path1W4 <= tff$W) & (path1W4 >= nw)))/(1004/100), digits = 1)

A2path1W4 <- round(sum(as.numeric((path1W4 >= tff$W)&(path1W4 >= nw)))/(1004/100), digits = 1)

A3path1W4 <- round(sum(as.numeric((path1W4 >= tff$W)&(path1W4 <= nw)))/(1004/100), digits=1)

A4path1W4 <- round(sum(as.numeric((path1W4 <= tff$W)&(path1W4 <= nw)))/(1004/100), digits = 1)


par(mfrow=c(2, 3))
plot(tff$W, path1W4, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 1',
     xlab= 'Path 1',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path1W4)
text(x=0.8, y=0.5, labels = A2path1W4)
text(x=0.5, y=0.2, labels = A3path1W4)
text(x=0.2, y=0.5, labels = A4path1W4)



A1path2W4 <- round(sum(as.numeric((path2W4 <= tff$W) & (path2W4 >= nw)))/(1004/100), digits = 1)

A2path2W4 <- round(sum(as.numeric((path2W4 >= tff$W)&(path2W4 >= nw)))/(1004/100), digits = 1)

A3path2W4 <- round(sum(as.numeric((path2W4 >= tff$W)&(path2W4 <= nw)))/(1004/100), digits=1)

A4path2W4 <- round(sum(as.numeric((path2W4 <= tff$W)&(path2W4 <= nw)))/(1004/100), digits = 1)

plot(tff$W, path2W4, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 2',
     xlab='Path 2',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path2W4)
text(x=0.8, y=0.5, labels = A2path2W4)
text(x=0.5, y=0.2, labels = A3path2W4)
text(x=0.2, y=0.5, labels = A4path2W4)

A1path3W4 <- round(sum(as.numeric((path3W4 <= tff$W) & (path3W4 >= nw)))/(1004/100), digits = 1)

A2path3W4 <- round(sum(as.numeric((path3W4 >= tff$W)&(path3W4 >= nw)))/(1004/100), digits = 1)

A3path3W4 <- round(sum(as.numeric((path3W4 >= tff$W)&(path3W4 <= nw)))/(1004/100), digits=1)

A4path3W4 <- round(sum(as.numeric((path3W4 <= tff$W)&(path3W4 <= nw)))/(1004/100), digits = 1)

plot(tff$W, path3W4, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main='Skewness path 3',
     xlab='Path 3',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path3W4)
text(x=0.8, y=0.5, labels = A2path3W4)
text(x=0.5, y=0.2, labels = A3path3W4)
text(x=0.2, y=0.5, labels = A4path3W4)

A1epsolW4 <- round(sum(as.numeric((epsolW4 <= tff$W) & (epsolW4 >= nw)))/(1004/100), digits = 1)

A2epsolW4 <- round(sum(as.numeric((epsolW4 >= tff$W)&(epsolW4 >= nw)))/(1004/100), digits = 1)

A3epsolW4 <- round(sum(as.numeric((epsolW4 >= tff$W)&(epsolW4 <= nw)))/(1004/100), digits=1)

A4epsolW4 <- round(sum(as.numeric((epsolW4 <= tff$W)&(epsolW4 <= nw)))/(1004/100), digits = 1)

isW4
plot(tff$W, epsolW4, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main=' Skewness eps',
     xlab='Solution term',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1epsolW4)
text(x=0.8, y=0.5, labels = A2epsolW4)
text(x=0.5, y=0.2, labels = A3epsolW4)
text(x=0.2, y=0.5, labels = A4epsolW4)

###### W5: frequency threshold 10, raw consistency threshold 0.9395######

ttW5 <- truthTable(tff, outcome="W", 
                   conditions = "SP, TP, OP, SM, CM", 
                   incl.cut=0.9395, n.cut=10, sort.by="incl, n", decreasing=TRUE, 
                   complete=FALSE, show.cases=FALSE)
ttW5


colnames(tff)

#parsimonious solution
psW5 <- eqmcc(ttW5, include="?", details=TRUE, show.cases=FALSE,  row.dom=TRUE, all.sol=FALSE)

psW5

#exclude contradictory assumptions & rows: tp*op + sp*op (contradict statement of necessity)

ttW5new <- ttW5
nec
rows <- findRows("tp*op + op*sp", ttW5new, remainders = FALSE)
rows

ttW5new$tt[as.character(rows), "OUT"] <- 0
ttW5new


#conservative solution
csW5 <- eqmcc(ttW5new, details=TRUE, PRI=TRUE, show.cases=FALSE, rowdom=TRUE, use.tilde=FALSE)
csW5

#enhanced parsimonious solution (all simplifying assumptions)

epsW5 <- eqmcc(ttW5new, include="?", details=TRUE, PRI=TRUE, show.cases=FALSE, 
               rowdom=TRUE, use.tilde=FALSE)
epsW5

#Z-Test for enhanced parsimonious solution, 0.8 (almost always sufficient)
# Calculate membership in solution term OP*CM + sp*TP*OP*SM
epsolW5 <- compute("OP*CM + sp*TP*OP*SM", data=tff)

obsepsolW5 <- as.numeric(epsolW5 <= tff$W)
zepsolW5 <- round(((((sum(obsepsolW5 == 1))/1004) - 0.8)-(1/(2*1004))) / (sqrt((0.8*(1-0.8))/1004)), digits=1)
zepsolW5
pzepsolW5 <- 2*pnorm(-abs(zepsolW5))
pzepsolW5


#intermediate solution with directional expectations SM -> W, CM -> W
isW5 <- eqmcc(ttW5new, details = TRUE, include = "?",  use.tilde = FALSE, PRI = TRUE,
              row.dom = TRUE, show.cases = FALSE, 
              dir.exp = "-, -, -, 1, 1")

isW5



#Test for skewness of sufficient conditions (enhanced parsimonious solution)
# calculate membership in paths of M1: OP*CM + sp*TP*OP*SM

path1W5 <- compute("OP*CM ", data=tff)
path2W5 <- compute("sp*TP*OP*SM ", data=tff)

A1path1W5 <- round(sum(as.numeric((path1W5 <= tff$W) & (path1W5 >= nw)))/(1004/100), digits = 1)

A2path1W5 <- round(sum(as.numeric((path1W5 >= tff$W)&(path1W5 >= nw)))/(1004/100), digits = 1)

A3path1W5 <- round(sum(as.numeric((path1W5 >= tff$W)&(path1W5 <= nw)))/(1004/100), digits=1)

A4path1W5 <- round(sum(as.numeric((path1W5 <= tff$W)&(path1W5 <= nw)))/(1004/100), digits = 1)


par(mfrow=c(2, 3))
plot(tff$W, path1W5, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 1',
     xlab= 'Path 1',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path1W5)
text(x=0.8, y=0.5, labels = A2path1W5)
text(x=0.5, y=0.2, labels = A3path1W5)
text(x=0.2, y=0.5, labels = A4path1W5)



A1path2W5 <- round(sum(as.numeric((path2W5 <= tff$W) & (path2W5 >= nw)))/(1004/100), digits = 1)

A2path2W5 <- round(sum(as.numeric((path2W5 >= tff$W)&(path2W5 >= nw)))/(1004/100), digits = 1)

A3path2W5 <- round(sum(as.numeric((path2W5 >= tff$W)&(path2W5 <= nw)))/(1004/100), digits=1)

A4path2W5 <- round(sum(as.numeric((path2W5 <= tff$W)&(path2W5 <= nw)))/(1004/100), digits = 1)

plot(tff$W, path2W5, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 2',
     xlab='Path 2',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path2W5)
text(x=0.8, y=0.5, labels = A2path2W5)
text(x=0.5, y=0.2, labels = A3path2W5)
text(x=0.2, y=0.5, labels = A4path2W5)

A1epsolW5 <- round(sum(as.numeric((epsolW5 <= tff$W) & (epsolW5 >= nw)))/(1004/100), digits = 1)

A2epsolW5 <- round(sum(as.numeric((epsolW5 >= tff$W)&(epsolW5 >= nw)))/(1004/100), digits = 1)

A3epsolW5 <- round(sum(as.numeric((epsolW5 >= tff$W)&(epsolW5 <= nw)))/(1004/100), digits=1)

A4epsolW5 <- round(sum(as.numeric((epsolW5 <= tff$W)&(epsolW5 <= nw)))/(1004/100), digits = 1)

isW5
plot(tff$W, epsolW5, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main=' Skewness eps',
     xlab='Solution term',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1epsolW5)
text(x=0.8, y=0.5, labels = A2epsolW5)
text(x=0.5, y=0.2, labels = A3epsolW5)
text(x=0.2, y=0.5, labels = A4epsolW5)



####W6 frequency threshold 10, raw consistency threshold 0.928####
ttW6 <- truthTable(tff, outcome="W", 
                   conditions = "SP, TP, OP, SM, CM", 
                   incl.cut=0.928, n.cut=10, sort.by="incl, n", decreasing=TRUE, 
                   complete=FALSE, show.cases=FALSE)
ttW6


colnames(tff)

#parsimonious solution
psW6 <- eqmcc(ttW6, include="?", details=TRUE, show.cases=FALSE,  row.dom=TRUE, all.sol=FALSE)

psW6

#exclude contradictory assumptions & rows: tp*op + sp*op (contradict statement of necessity)

ttW6new <- ttW6
nec
rows <- findRows("tp*op + op*sp", ttW6new, remainders = FALSE)
rows

ttW6new$tt[as.character(rows), "OUT"] <- 0
ttW6new


#conservative solution
csW6 <- eqmcc(ttW6new, details=TRUE, PRI=TRUE, show.cases=FALSE, rowdom=TRUE, use.tilde=FALSE)
csW6

#enhanced parsimonious solution (all simplifying assumptions)

epsW6 <- eqmcc(ttW6new, include="?", details=TRUE, PRI=TRUE, show.cases=FALSE, 
               rowdom=TRUE, use.tilde=FALSE)
epsW6


#Z-Test for enhanced parsimonious solution, 0.8 (almost always sufficient)
# Calculate membership in solution term OP*CM + OP*SM
epsolW6 <- compute("OP*CM + OP*SM", data=tff)

obsepsolW6 <- as.numeric(epsolW6 <= tff$W)
zepsolW6 <- round(((((sum(obsepsolW6 == 1))/1004) - 0.8)-(1/(2*1004))) / (sqrt((0.8*(1-0.8))/1004)), digits=1)
zepsolW6
pzepsolW6 <- 2*pnorm(-abs(zepsolW6))
pzepsolW6


#intermediate solution with directional expectations SM -> W, CM -> W
isW6 <- eqmcc(ttW6new, details = TRUE, include = "?",  use.tilde = FALSE, PRI = TRUE,
              row.dom = TRUE, show.cases = FALSE, 
              dir.exp = "-, -, -, 1, 1")

isW6

#Test for skewness of sufficient conditions (enhanced parsimonious solution)
# calculate membership in paths of M1: OP*CM + OP*SM

path1W6 <- compute("OP*CM ", data=tff)
path2W6 <- compute("OP*SM ", data=tff)

A1path1W6 <- round(sum(as.numeric((path1W6 <= tff$W) & (path1W6 >= nw)))/(1004/100), digits = 1)

A2path1W6 <- round(sum(as.numeric((path1W6 >= tff$W)&(path1W6 >= nw)))/(1004/100), digits = 1)

A3path1W6 <- round(sum(as.numeric((path1W6 >= tff$W)&(path1W6 <= nw)))/(1004/100), digits=1)

A4path1W6 <- round(sum(as.numeric((path1W6 <= tff$W)&(path1W6 <= nw)))/(1004/100), digits = 1)


par(mfrow=c(2, 3))
plot(tff$W, path1W6, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 1',
     xlab= 'Path 1',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path1W6)
text(x=0.8, y=0.5, labels = A2path1W6)
text(x=0.5, y=0.2, labels = A3path1W6)
text(x=0.2, y=0.5, labels = A4path1W6)



A1path2W6 <- round(sum(as.numeric((path2W6 <= tff$W) & (path2W6 >= nw)))/(1004/100), digits = 1)

A2path2W6 <- round(sum(as.numeric((path2W6 >= tff$W)&(path2W6 >= nw)))/(1004/100), digits = 1)

A3path2W6 <- round(sum(as.numeric((path2W6 >= tff$W)&(path2W6 <= nw)))/(1004/100), digits=1)

A4path2W6 <- round(sum(as.numeric((path2W6 <= tff$W)&(path2W6 <= nw)))/(1004/100), digits = 1)

plot(tff$W, path2W6, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 2',
     xlab='Path 2',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path2W6)
text(x=0.8, y=0.5, labels = A2path2W6)
text(x=0.5, y=0.2, labels = A3path2W6)
text(x=0.2, y=0.5, labels = A4path2W6)

A1epsolW6 <- round(sum(as.numeric((epsolW6 <= tff$W) & (epsolW6 >= nw)))/(1004/100), digits = 1)

A2epsolW6 <- round(sum(as.numeric((epsolW6 >= tff$W)&(epsolW6 >= nw)))/(1004/100), digits = 1)

A3epsolW6 <- round(sum(as.numeric((epsolW6 >= tff$W)&(epsolW6 <= nw)))/(1004/100), digits=1)

A4epsolW6 <- round(sum(as.numeric((epsolW6 <= tff$W)&(epsolW6 <= nw)))/(1004/100), digits = 1)

isW6
plot(tff$W, epsolW6, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main=' Skewness eps',
     xlab='Solution term',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1epsolW6)
text(x=0.8, y=0.5, labels = A2epsolW6)
text(x=0.5, y=0.2, labels = A3epsolW6)
text(x=0.2, y=0.5, labels = A4epsolW6)


#####W7 frequency threshold 50, raw consistency threshold 0.928####

ttW7 <- truthTable(tff, outcome="W", 
                   conditions = "SP, TP, OP, SM, CM", 
                   incl.cut=0.928, n.cut=50, sort.by="incl, n", decreasing=TRUE, 
                   complete=FALSE, show.cases=FALSE)
ttW7


colnames(tff)

#parsimonious solution
psW7 <- eqmcc(ttW7, include="?", details=TRUE, show.cases=FALSE,  row.dom=TRUE, all.sol=FALSE)

psW7

#exclude contradictory assumptions & rows: tp*op + sp*op (contradict statement of necessity)

ttW7new <- ttW7
nec
rows <- findRows("tp*op + op*sp", ttW7new, remainders = FALSE)
rows

ttW7new$tt[as.character(rows), "OUT"] <- 0
ttW7new


#conservative solution
csW7 <- eqmcc(ttW7new, details=TRUE, PRI=TRUE, show.cases=FALSE, rowdom=TRUE, use.tilde=FALSE)
csW7

#enhanced parsimonious solution (all simplifying assumptions)

epsW7 <- eqmcc(ttW7new, include="?", details=TRUE, PRI=TRUE, show.cases=FALSE, 
               rowdom=TRUE, use.tilde=FALSE)
epsW7


#Z-Test for enhanced parsimonious solution, 0.8 (almost always sufficient)
# Calculate membership in solution term SP*TP*SM
epsolW7 <- compute("SP*TP*SM", data=tff)

obsepsolW7 <- as.numeric(epsolW7 <= tff$W)
zepsolW7 <- round(((((sum(obsepsolW7 == 1))/1004) - 0.8)-(1/(2*1004))) / (sqrt((0.8*(1-0.8))/1004)), digits=1)
zepsolW7
pzepsolW7 <- 2*pnorm(-abs(zepsolW7))
pzepsolW7


#intermediate solution with directional expectations SM -> W, CM -> W
isW7 <- eqmcc(ttW7new, details = TRUE, include = "?",  use.tilde = FALSE, PRI = TRUE,
              row.dom = TRUE, show.cases = FALSE, 
              dir.exp = "-, -, -, 1, 1")

isW7

#Test for skewness of sufficient conditions (enhanced parsimonious solution)
# calculate membership in paths of M1: SP*TP*SM

path1W7 <- compute("SP*TP*SM", data=tff)

A1path1W7 <- round(sum(as.numeric((path1W7 <= tff$W) & (path1W7 >= nw)))/(1004/100), digits = 1)

A2path1W7 <- round(sum(as.numeric((path1W7 >= tff$W)&(path1W7 >= nw)))/(1004/100), digits = 1)

A3path1W7 <- round(sum(as.numeric((path1W7 >= tff$W)&(path1W7 <= nw)))/(1004/100), digits=1)

A4path1W7 <- round(sum(as.numeric((path1W7 <= tff$W)&(path1W7 <= nw)))/(1004/100), digits = 1)


par(mfrow=c(2, 3))
plot(tff$W, path1W7, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 1',
     xlab= 'Path 1',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path1W7)
text(x=0.8, y=0.5, labels = A2path1W7)
text(x=0.5, y=0.2, labels = A3path1W7)
text(x=0.2, y=0.5, labels = A4path1W7)


A1epsolW7 <- round(sum(as.numeric((epsolW7 <= tff$W) & (epsolW7 >= nw)))/(1004/100), digits = 1)

A2epsolW7 <- round(sum(as.numeric((epsolW7 >= tff$W)&(epsolW7 >= nw)))/(1004/100), digits = 1)

A3epsolW7 <- round(sum(as.numeric((epsolW7 >= tff$W)&(epsolW7 <= nw)))/(1004/100), digits=1)

A4epsolW7 <- round(sum(as.numeric((epsolW7 <= tff$W)&(epsolW7 <= nw)))/(1004/100), digits = 1)

isW7
plot(tff$W, epsolW7, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main=' Skewness eps',
     xlab='Solution term',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1epsolW7)
text(x=0.8, y=0.5, labels = A2epsolW7)
text(x=0.5, y=0.2, labels = A3epsolW7)
text(x=0.2, y=0.5, labels = A4epsolW7)

#####W8 frequency threshold 50, raw consistency threshold 0.858####

ttW8 <- truthTable(tff, outcome="W", 
                   conditions = "SP, TP, OP, SM, CM", 
                   incl.cut=0.858, n.cut=50, sort.by="incl, n", decreasing=TRUE, 
                   complete=FALSE, show.cases=FALSE)
ttW8


colnames(tff)

#parsimonious solution
psW8 <- eqmcc(ttW8, include="?", details=TRUE, show.cases=FALSE,  row.dom=TRUE, all.sol=FALSE)

psW8

#exclude contradictory assumptions & rows: tp*op + sp*op (contradict statement of necessity)

ttW8new <- ttW8
nec
rows <- findRows("tp*op + op*sp", ttW8new, remainders = FALSE)
rows

ttW8new$tt[as.character(rows), "OUT"] <- 0
ttW8new


#conservative solution
csW8 <- eqmcc(ttW8new, details=TRUE, PRI=TRUE, show.cases=FALSE, rowdom=TRUE, use.tilde=FALSE)
csW8

#enhanced parsimonious solution (all simplifying assumptions)

epsW8 <- eqmcc(ttW8new, include="?", details=TRUE, PRI=TRUE, show.cases=FALSE, 
               rowdom=TRUE, use.tilde=FALSE)
epsW8


#Z-Test for enhanced parsimonious solution, 0.8 (almost always sufficient)
# Calculate membership in solution term SP*OP
epsolW8 <- compute("SP*OP", data=tff)

obsepsolW8 <- as.numeric(epsolW8 <= tff$W)
zepsolW8 <- round(((((sum(obsepsolW8 == 1))/1004) - 0.8)-(1/(2*1004))) / (sqrt((0.8*(1-0.8))/1004)), digits=1)
zepsolW8
pzepsolW8 <- 2*pnorm(-abs(zepsolW8))
pzepsolW8


#intermediate solution with directional expectations SM -> W, CM -> W
isW8 <- eqmcc(ttW8new, details = TRUE, include = "?",  use.tilde = FALSE, PRI = TRUE,
              row.dom = TRUE, show.cases = FALSE, 
              dir.exp = "-, -, -, 1, 1")

isW8

#Test for skewness of sufficient conditions (enhanced parsimonious solution)
# calculate membership in paths of M1: SP*OP

path1W8 <- compute("SP*OP ", data=tff)

A1path1W8 <- round(sum(as.numeric((path1W8 <= tff$W) & (path1W8 >= nw)))/(1004/100), digits = 1)

A2path1W8 <- round(sum(as.numeric((path1W8 >= tff$W)&(path1W8 >= nw)))/(1004/100), digits = 1)

A3path1W8 <- round(sum(as.numeric((path1W8 >= tff$W)&(path1W8 <= nw)))/(1004/100), digits=1)

A4path1W8 <- round(sum(as.numeric((path1W8 <= tff$W)&(path1W8 <= nw)))/(1004/100), digits = 1)


par(mfrow=c(2, 3))
plot(tff$W, path1W8, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 1',
     xlab= 'Path 1',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path1W8)
text(x=0.8, y=0.5, labels = A2path1W8)
text(x=0.5, y=0.2, labels = A3path1W8)
text(x=0.2, y=0.5, labels = A4path1W8)


A1epsolW8 <- round(sum(as.numeric((epsolW8 <= tff$W) & (epsolW8 >= nw)))/(1004/100), digits = 1)

A2epsolW8 <- round(sum(as.numeric((epsolW8 >= tff$W)&(epsolW8 >= nw)))/(1004/100), digits = 1)

A3epsolW8 <- round(sum(as.numeric((epsolW8 >= tff$W)&(epsolW8 <= nw)))/(1004/100), digits=1)

A4epsolW8 <- round(sum(as.numeric((epsolW8 <= tff$W)&(epsolW8 <= nw)))/(1004/100), digits = 1)

isW8
plot(tff$W, epsolW8, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main=' Skewness eps',
     xlab='Solution term',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1epsolW8)
text(x=0.8, y=0.5, labels = A2epsolW8)
text(x=0.5, y=0.2, labels = A3epsolW8)
text(x=0.2, y=0.5, labels = A4epsolW8)

#####W9 frequency threshold 50, raw consistency threshold 0.839####

ttW9 <- truthTable(tff, outcome="W", 
                   conditions = "SP, TP, OP, SM, CM", 
                   incl.cut=0.839, n.cut=50, sort.by="incl, n", decreasing=TRUE, 
                   complete=FALSE, show.cases=FALSE)
ttW9


colnames(tff)

#parsimonious solution
psW9 <- eqmcc(ttW9, include="?", details=TRUE, show.cases=FALSE,  row.dom=TRUE, all.sol=FALSE)

psW9

#exclude contradictory assumptions & rows: tp*op + sp*op (contradict statement of necessity)

ttW9new <- ttW9
nec
rows <- findRows("tp*op + op*sp", ttW9new, remainders = FALSE)
rows

ttW9new$tt[as.character(rows), "OUT"] <- 0
ttW9new


#conservative solution
csW9 <- eqmcc(ttW9new, details=TRUE, PRI=TRUE, show.cases=FALSE, rowdom=TRUE, use.tilde=FALSE)
csW9

#enhanced parsimonious solution (all simplifying assumptions)

epsW9 <- eqmcc(ttW9new, include="?", details=TRUE, PRI=TRUE, show.cases=FALSE, 
               rowdom=TRUE, use.tilde=FALSE)
epsW9

#Z-Test for enhanced parsimonious solution, 0.8 (almost always sufficient)
# Calculate membership in solution term SP*OP + TP*OP
epsolW9 <- compute("SP*OP + TP*OP", data=tff)

obsepsolW9 <- as.numeric(epsolW9 <= tff$W)
zepsolW9 <- round(((((sum(obsepsolW9 == 1))/1004) - 0.8)-(1/(2*1004))) / (sqrt((0.8*(1-0.8))/1004)), digits=1)
zepsolW9
pzepsolW9 <- 2*pnorm(-abs(zepsolW9))
pzepsolW9


#intermediate solution with directional expectations SM -> W, CM -> W
isW9 <- eqmcc(ttW9new, details = TRUE, include = "?",  use.tilde = FALSE, PRI = TRUE,
              row.dom = TRUE, show.cases = FALSE, 
              dir.exp = "-, -, -, 1, 1")

isW9

#Test for skewness of sufficient conditions (enhanced parsimonious solution)
# calculate membership in paths of M1: SP*OP + TP*OP

path1W9 <- compute("SP*OP ", data=tff)

A1path1W9 <- round(sum(as.numeric((path1W9 <= tff$W) & (path1W9 >= nw)))/(1004/100), digits = 1)

A2path1W9 <- round(sum(as.numeric((path1W9 >= tff$W)&(path1W9 >= nw)))/(1004/100), digits = 1)

A3path1W9 <- round(sum(as.numeric((path1W9 >= tff$W)&(path1W9 <= nw)))/(1004/100), digits=1)

A4path1W9 <- round(sum(as.numeric((path1W9 <= tff$W)&(path1W9 <= nw)))/(1004/100), digits = 1)


par(mfrow=c(2, 3))
plot(tff$W, path1W9, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 1',
     xlab= 'Path 1',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path1W9)
text(x=0.8, y=0.5, labels = A2path1W9)
text(x=0.5, y=0.2, labels = A3path1W9)
text(x=0.2, y=0.5, labels = A4path1W9)

path2W9 <- compute("TP*OP ", data=tff)

A1path2W9 <- round(sum(as.numeric((path2W9 <= tff$W) & (path2W9 >= nw)))/(1004/100), digits = 1)

A2path2W9 <- round(sum(as.numeric((path2W9 >= tff$W)&(path2W9 >= nw)))/(1004/100), digits = 1)

A3path2W9 <- round(sum(as.numeric((path2W9 >= tff$W)&(path2W9 <= nw)))/(1004/100), digits=1)

A4path2W9 <- round(sum(as.numeric((path2W9 <= tff$W)&(path2W9 <= nw)))/(1004/100), digits = 1)


plot(tff$W, path2W9, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 2',
     xlab= 'Path 2',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path2W9)
text(x=0.8, y=0.5, labels = A2path2W9)
text(x=0.5, y=0.2, labels = A3path2W9)
text(x=0.2, y=0.5, labels = A4path2W9)



A1epsolW9 <- round(sum(as.numeric((epsolW9 <= tff$W) & (epsolW9 >= nw)))/(1004/100), digits = 1)

A2epsolW9 <- round(sum(as.numeric((epsolW9 >= tff$W)&(epsolW9 >= nw)))/(1004/100), digits = 1)

A3epsolW9 <- round(sum(as.numeric((epsolW9 >= tff$W)&(epsolW9 <= nw)))/(1004/100), digits=1)

A4epsolW9 <- round(sum(as.numeric((epsolW9 <= tff$W)&(epsolW9 <= nw)))/(1004/100), digits = 1)

isW9
plot(tff$W, epsolW9, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main=' Skewness eps',
     xlab='Solution term',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1epsolW9)
text(x=0.8, y=0.5, labels = A2epsolW9)
text(x=0.5, y=0.2, labels = A3epsolW9)
text(x=0.2, y=0.5, labels = A4epsolW9)


#######################ANALYSIS OF SUFFICIENCY FOR w ####################


tff <- read.csv("tff_d1_c2.csv", row.names=1, header = TRUE, sep = ";",  dec = ",")
describe(tff)

nw <- 1-tff$W

# unenable assumptions: 

suf <- "OP*CM + SP*TP*CM + sp*TP*OP*SM + SP*TP*op*SM"
suf


#####w1  frequency threshold 5, raw consistency threshold 0.945 #####


ttw1 <- truthTable(tff, outcome="~W", 
                   conditions = "SP, TP, OP, SM, CM", 
                   incl.cut=0.945, n.cut=5, sort.by="incl, n", decreasing=TRUE, 
                   complete=FALSE, show.cases=FALSE)
ttw1


#parsimonious solution
psw1 <- eqmcc(ttw1, include="?", details=TRUE, show.cases=FALSE,  row.dom=TRUE, all.sol=FALSE)

psw1

#exclude contradictory assumptions & rows: OP*CM + SP*TP*CM + sp*TP*OP*SM + SP*TP*op*SM (contradict statement of sufficiency)

ttw1new <- ttw1
suf
rows <- findRows("OP*CM + SP*TP*CM + sp*TP*OP*SM + SP*TP*op*SM", ttw1new, remainders = FALSE)
rows

ttw1new$tt[as.character(rows), "OUT"] <- 0
ttw1new



#conservative solution
csw1 <- eqmcc(ttw1new, details=TRUE, PRI=TRUE, show.cases=FALSE, rowdom=TRUE, use.tilde=FALSE)
csw1

#enhanced parsimonious solution (all simplifying assumptions)

epsw1 <- eqmcc(ttw1new, include="?", details=TRUE, PRI=TRUE, show.cases=FALSE, 
               rowdom=TRUE, use.tilde=FALSE)
epsw1


#Z-Test for enhanced parsimonious solution, 0.8 (almost always sufficient)
# Calculate membership in solution term  tp*op*CM + sp*op*sm*CM
epsolw1 <- compute("tp*op*CM + sp*op*sm*CM ", data=tff)

obsepsolw1 <- as.numeric(epsolw1 <= nw)
zepsolw1 <- round(((((sum(obsepsolw1 == 1))/1004) - 0.8)-(1/(2*1004))) / (sqrt((0.8*(1-0.8))/1004)), digits=1)
zepsolw1
pzepsolw1 <- 2*pnorm(-abs(zepsolw1))
pzepsolw1



#intermediate solution with Directional expectations: sp->w, tp->w, op->w.

isw1 <- eqmcc(ttw1new, details = TRUE, include = "?",  use.tilde = FALSE, PRI = TRUE,
              row.dom = TRUE, show.cases = FALSE, 
              dir.exp = "0, 0, 0, -, -")

isw1


#Test for skewness of sufficient conditions (enhanced parsimonious solution)
# calculate membership in paths of M1: tp*op*CM + sp*op*sm*CM
path1w1 <- compute("tp*op*CM ", data=tff)
path2w1 <- compute("sp*op*sm*CM ", data=tff)


A1path1w1 <- round(sum(as.numeric((path1w1 <= nw) & (path1w1 >= tff$W)))/(1004/100), digits = 1)

A2path1w1 <- round(sum(as.numeric((path1w1 >= nw)&(path1w1 >= tff$W)))/(1004/100), digits = 1)

A3path1w1 <- round(sum(as.numeric((path1w1 >= nw)&(path1w1 <= tff$W)))/(1004/100), digits=1)

A4path1w1 <- round(sum(as.numeric((path1w1 <= nw)&(path1w1 <= tff$W)))/(1004/100), digits = 1)


par(mfrow=c(2, 3))
plot(nw, path1w1, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 1',
     xlab= 'Path 1',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path1w1)
text(x=0.8, y=0.5, labels = A2path1w1)
text(x=0.5, y=0.2, labels = A3path1w1)
text(x=0.2, y=0.5, labels = A4path1w1)



A1path2w1 <- round(sum(as.numeric((path2w1 <= nw) & (path2w1 >= tff$W)))/(1004/100), digits = 1)

A2path2w1 <- round(sum(as.numeric((path2w1 >= nw)&(path2w1 >= tff$W)))/(1004/100), digits = 1)

A3path2w1 <- round(sum(as.numeric((path2w1 >= nw)&(path2w1 <= tff$W)))/(1004/100), digits=1)

A4path2w1 <- round(sum(as.numeric((path2w1 <= nw)&(path2w1 <= tff$W)))/(1004/100), digits = 1)


plot(nw, path2w1, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 2',
     xlab= 'Path 2',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path2w1)
text(x=0.8, y=0.5, labels = A2path2w1)
text(x=0.5, y=0.2, labels = A3path2w1)
text(x=0.2, y=0.5, labels = A4path2w1)


A1epsolw1 <- round(sum(as.numeric((epsolw1 <= nw) & (epsolw1 >= tff$W)))/(1004/100), digits = 1)

A2epsolw1 <- round(sum(as.numeric((epsolw1 >= nw)&(epsolw1 >= tff$W)))/(1004/100), digits = 1)

A3epsolw1 <- round(sum(as.numeric((epsolw1 >= nw)&(epsolw1 <= tff$W)))/(1004/100), digits=1)

A4epsolw1 <- round(sum(as.numeric((epsolw1 <= nw)&(epsolw1 <= tff$W)))/(1004/100), digits = 1)


plot(nw, epsolw1, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness eps',
     xlab= 'solution term',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1epsolw1)
text(x=0.8, y=0.5, labels = A2epsolw1)
text(x=0.5, y=0.2, labels = A3epsolw1)
text(x=0.2, y=0.5, labels = A4epsolw1)



#####w2  frequency threshold 5, raw consistency threshold 0.943 #####


ttw2 <- truthTable(tff, outcome="~W", 
                   conditions = "SP, TP, OP, SM, CM", 
                   incl.cut=0.943, n.cut=5, sort.by="incl, n", decreasing=TRUE, 
                   complete=FALSE, show.cases=FALSE)
ttw2

#parsimonious solution
psw2 <- eqmcc(ttw2, include="?", details=TRUE, show.cases=FALSE,  row.dom=TRUE, all.sol=FALSE)

psw2

#exclude contradictory assumptions & rows: OP*CM + SP*TP*CM + sp*TP*OP*SM + SP*TP*op*SM (contradict statement of sufficiency)

ttw2new <- ttw2
suf
rows <- findRows("OP*CM + SP*TP*CM + sp*TP*OP*SM + SP*TP*op*SM", ttw2new, remainders = FALSE)
rows

ttw2new$tt[as.character(rows), "OUT"] <- 0
ttw2new


#conservative solution
csw2 <- eqmcc(ttw2new, details=TRUE, PRI=TRUE, show.cases=FALSE, rowdom=TRUE, use.tilde=FALSE)
csw2

#enhanced parsimonious solution (all simplifying assumptions)

epsw2 <- eqmcc(ttw2new, include="?", details=TRUE, PRI=TRUE, show.cases=FALSE, 
               rowdom=TRUE, use.tilde=FALSE)
epsw2

#Z-Test for enhanced parsimonious solution, 0.8 (almost always sufficient)
# Calculate membership in solution term  tp*op*CM + sp*op*sm*CM + sp*TP*op*SM*cm
epsolw2 <- compute("tp*op*CM + sp*op*sm*CM + sp*TP*op*SM*cm ", data=tff)

obsepsolw2 <- as.numeric(epsolw2 <= nw)
zepsolw2 <- round(((((sum(obsepsolw2 == 1))/1004) - 0.8)-(1/(2*1004))) / (sqrt((0.8*(1-0.8))/1004)), digits=1)
zepsolw2
pzepsolw2 <- 2*pnorm(-abs(zepsolw2))
pzepsolw2

#intermediate solution with Directional expectations: sp->w, tp->w, op->w.
isw2 <- eqmcc(ttw2new, details = TRUE, include = "?",  use.tilde = FALSE, PRI = TRUE,
              row.dom = TRUE, show.cases = FALSE, 
              dir.exp = "0, 0, 0, -, -")

isw2


#Test for skewness of sufficient conditions (enhanced parsimonious solution)
# calculate membership in paths of M1: tp*op*CM + sp*op*sm*CM + sp*TP*op*SM*cm
path1w2 <- compute("tp*op*CM ", data=tff)
path2w2 <- compute("sp*op*sm*CM ", data=tff)
path3w2 <- compute("sp*TP*op*SM*cm ", data=tff)

A1path1w2 <- round(sum(as.numeric((path1w2 <= nw) & (path1w2 >= tff$W)))/(1004/100), digits = 1)

A2path1w2 <- round(sum(as.numeric((path1w2 >= nw)&(path1w2 >= tff$W)))/(1004/100), digits = 1)

A3path1w2 <- round(sum(as.numeric((path1w2 >= nw)&(path1w2 <= tff$W)))/(1004/100), digits=1)

A4path1w2 <- round(sum(as.numeric((path1w2 <= nw)&(path1w2 <= tff$W)))/(1004/100), digits = 1)


par(mfrow=c(2, 3))
plot(nw, path1w2, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 1',
     xlab= 'Path 1',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path1w2)
text(x=0.8, y=0.5, labels = A2path1w2)
text(x=0.5, y=0.2, labels = A3path1w2)
text(x=0.2, y=0.5, labels = A4path1w2)



A1path2w2 <- round(sum(as.numeric((path2w2 <= nw) & (path2w2 >= tff$W)))/(1004/100), digits = 1)

A2path2w2 <- round(sum(as.numeric((path2w2 >= nw)&(path2w2 >= tff$W)))/(1004/100), digits = 1)

A3path2w2 <- round(sum(as.numeric((path2w2 >= nw)&(path2w2 <= tff$W)))/(1004/100), digits=1)

A4path2w2 <- round(sum(as.numeric((path2w2 <= nw)&(path2w2 <= tff$W)))/(1004/100), digits = 1)


plot(nw, path2w2, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 2',
     xlab= 'Path 2',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path2w2)
text(x=0.8, y=0.5, labels = A2path2w2)
text(x=0.5, y=0.2, labels = A3path2w2)
text(x=0.2, y=0.5, labels = A4path2w2)

A1path3w2 <- round(sum(as.numeric((path3w2 <= nw) & (path3w2 >= tff$W)))/(1004/100), digits = 1)

A2path3w2 <- round(sum(as.numeric((path3w2 >= nw)&(path3w2 >= tff$W)))/(1004/100), digits = 1)

A3path3w2 <- round(sum(as.numeric((path3w2 >= nw)&(path3w2 <= tff$W)))/(1004/100), digits=1)

A4path3w2 <- round(sum(as.numeric((path3w2 <= nw)&(path3w2 <= tff$W)))/(1004/100), digits = 1)

plot(nw, path3w2, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 3',
     xlab= 'Path 3',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path3w2)
text(x=0.8, y=0.5, labels = A2path3w2)
text(x=0.5, y=0.2, labels = A3path3w2)
text(x=0.2, y=0.5, labels = A4path3w2)


A1epsolw2 <- round(sum(as.numeric((epsolw2 <= nw) & (epsolw2 >= tff$W)))/(1004/100), digits = 1)

A2epsolw2 <- round(sum(as.numeric((epsolw2 >= nw)&(epsolw2 >= tff$W)))/(1004/100), digits = 1)

A3epsolw2 <- round(sum(as.numeric((epsolw2 >= nw)&(epsolw2 <= tff$W)))/(1004/100), digits=1)

A4epsolw2 <- round(sum(as.numeric((epsolw2 <= nw)&(epsolw2 <= tff$W)))/(1004/100), digits = 1)


plot(nw, epsolw2, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness eps',
     xlab= 'solution term',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1epsolw2)
text(x=0.8, y=0.5, labels = A2epsolw2)
text(x=0.5, y=0.2, labels = A3epsolw2)
text(x=0.2, y=0.5, labels = A4epsolw2)

#####w3  frequency threshold 5, raw consistency threshold 0.939 #####


ttw3 <- truthTable(tff, outcome="~W", 
                   conditions = "SP, TP, OP, SM, CM", 
                   incl.cut=0.939, n.cut=5, sort.by="incl, n", decreasing=TRUE, 
                   complete=FALSE, show.cases=FALSE)
ttw3

#parsimonious solution
psw3 <- eqmcc(ttw3, include="?", details=TRUE, show.cases=FALSE,  row.dom=TRUE, all.sol=FALSE)

psw3

#exclude contradictory assumptions & rows: OP*CM + SP*TP*CM + sp*TP*OP*SM + SP*TP*op*SM (contradict statement of sufficiency)

ttw3new <- ttw3
suf
rows <- findRows("OP*CM + SP*TP*CM + sp*TP*OP*SM + SP*TP*op*SM", ttw3new, remainders = FALSE)
rows

ttw3new$tt[as.character(rows), "OUT"] <- 0
ttw3new


#conservative solution
csw3 <- eqmcc(ttw3new, details=TRUE, PRI=TRUE, show.cases=FALSE, rowdom=TRUE, use.tilde=FALSE)
csw3

#enhanced parsimonious solution (all simplifying assumptions)

epsw3 <- eqmcc(ttw3new, include="?", details=TRUE, PRI=TRUE, show.cases=FALSE, 
               rowdom=TRUE, use.tilde=FALSE)
epsw3


#Z-Test for enhanced parsimonious solution, 0.8 (almost always sufficient)
# Calculate membership in solution term  sp*op*CM + sp*op*SM + tp*op*CM + tp*op*SM
epsolw3 <- compute("sp*op*CM + sp*op*SM + tp*op*CM + tp*op*SM ", data=tff)

obsepsolw3 <- as.numeric(epsolw3 <= nw)
zepsolw3 <- round(((((sum(obsepsolw3 == 1))/1004) - 0.8)-(1/(2*1004))) / (sqrt((0.8*(1-0.8))/1004)), digits=1)
zepsolw3
pzepsolw3 <- 2*pnorm(-abs(zepsolw3))
pzepsolw3


#intermediate solution with Directional expectations: sp->w, tp->w, op->w.
isw3 <- eqmcc(ttw3new, details = TRUE, include = "?",  use.tilde = FALSE, PRI = TRUE,
              row.dom = TRUE, show.cases = FALSE, 
              dir.exp = "0, 0, 0, -, -")

isw3

#Test for skewness of sufficient conditions (enhanced parsimonious solution)
# calculate membership in paths of M1: sp*op*CM + sp*op*SM + tp*op*CM + tp*op*SM
path1w3 <- compute("sp*op*CM ", data=tff)
path2w3 <- compute("sp*op*SM ", data=tff)
path3w3 <- compute("tp*op*CM ", data=tff)
path4w3 <- compute("tp*op*SM ", data=tff)

A1path1w3 <- round(sum(as.numeric((path1w3 <= nw) & (path1w3 >= tff$W)))/(1004/100), digits = 1)

A2path1w3 <- round(sum(as.numeric((path1w3 >= nw)&(path1w3 >= tff$W)))/(1004/100), digits = 1)

A3path1w3 <- round(sum(as.numeric((path1w3 >= nw)&(path1w3 <= tff$W)))/(1004/100), digits=1)

A4path1w3 <- round(sum(as.numeric((path1w3 <= nw)&(path1w3 <= tff$W)))/(1004/100), digits = 1)


par(mfrow=c(2, 3))
plot(nw, path1w3, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 1',
     xlab= 'Path 1',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path1w3)
text(x=0.8, y=0.5, labels = A2path1w3)
text(x=0.5, y=0.2, labels = A3path1w3)
text(x=0.2, y=0.5, labels = A4path1w3)



A1path2w3 <- round(sum(as.numeric((path2w3 <= nw) & (path2w3 >= tff$W)))/(1004/100), digits = 1)

A2path2w3 <- round(sum(as.numeric((path2w3 >= nw)&(path2w3 >= tff$W)))/(1004/100), digits = 1)

A3path2w3 <- round(sum(as.numeric((path2w3 >= nw)&(path2w3 <= tff$W)))/(1004/100), digits=1)

A4path2w3 <- round(sum(as.numeric((path2w3 <= nw)&(path2w3 <= tff$W)))/(1004/100), digits = 1)


plot(nw, path2w3, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 2',
     xlab= 'Path 2',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path2w3)
text(x=0.8, y=0.5, labels = A2path2w3)
text(x=0.5, y=0.2, labels = A3path2w3)
text(x=0.2, y=0.5, labels = A4path2w3)

A1path3w3 <- round(sum(as.numeric((path3w3 <= nw) & (path3w3 >= tff$W)))/(1004/100), digits = 1)

A2path3w3 <- round(sum(as.numeric((path3w3 >= nw)&(path3w3 >= tff$W)))/(1004/100), digits = 1)

A3path3w3 <- round(sum(as.numeric((path3w3 >= nw)&(path3w3 <= tff$W)))/(1004/100), digits=1)

A4path3w3 <- round(sum(as.numeric((path3w3 <= nw)&(path3w3 <= tff$W)))/(1004/100), digits = 1)

plot(nw, path3w3, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 3',
     xlab= 'Path 3',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path3w3)
text(x=0.8, y=0.5, labels = A2path3w3)
text(x=0.5, y=0.2, labels = A3path3w3)
text(x=0.2, y=0.5, labels = A4path3w3)

A1path4w3 <- round(sum(as.numeric((path4w3 <= nw) & (path4w3 >= tff$W)))/(1004/100), digits = 1)

A2path4w3 <- round(sum(as.numeric((path4w3 >= nw)&(path4w3 >= tff$W)))/(1004/100), digits = 1)

A3path4w3 <- round(sum(as.numeric((path4w3 >= nw)&(path4w3 <= tff$W)))/(1004/100), digits=1)

A4path4w3 <- round(sum(as.numeric((path4w3 <= nw)&(path4w3 <= tff$W)))/(1004/100), digits = 1)


plot(nw, path4w3, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 4',
     xlab= 'Path 4',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path4w3)
text(x=0.8, y=0.5, labels = A2path4w3)
text(x=0.5, y=0.2, labels = A3path4w3)
text(x=0.2, y=0.5, labels = A4path4w3)



A1epsolw3 <- round(sum(as.numeric((epsolw3 <= nw) & (epsolw3 >= tff$W)))/(1004/100), digits = 1)

A2epsolw3 <- round(sum(as.numeric((epsolw3 >= nw)&(epsolw3 >= tff$W)))/(1004/100), digits = 1)

A3epsolw3 <- round(sum(as.numeric((epsolw3 >= nw)&(epsolw3 <= tff$W)))/(1004/100), digits=1)

A4epsolw3 <- round(sum(as.numeric((epsolw3 <= nw)&(epsolw3 <= tff$W)))/(1004/100), digits = 1)


plot(nw, epsolw3, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness eps',
     xlab= 'solution term',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1epsolw3)
text(x=0.8, y=0.5, labels = A2epsolw3)
text(x=0.5, y=0.2, labels = A3epsolw3)
text(x=0.2, y=0.5, labels = A4epsolw3)

#####w4  frequency threshold 10, raw consistency threshold 0.945 #####


ttw4 <- truthTable(tff, outcome="~W", 
                   conditions = "SP, TP, OP, SM, CM", 
                   incl.cut=0.945, n.cut=10, sort.by="incl, n", decreasing=TRUE, 
                   complete=FALSE, show.cases=FALSE)
ttw4

#parsimonious solution
psw4 <- eqmcc(ttw4, include="?", details=TRUE, show.cases=FALSE,  row.dom=TRUE, all.sol=FALSE)

psw4

#exclude contradictory assumptions & rows: OP*CM + SP*TP*CM + sp*TP*OP*SM + SP*TP*op*SM (contradict statement of sufficiency)

ttw4new <- ttw4
suf
rows <- findRows("OP*CM + SP*TP*CM + sp*TP*OP*SM + SP*TP*op*SM", ttw4new, remainders = FALSE)
rows

ttw4new$tt[as.character(rows), "OUT"] <- 0
ttw4new


#conservative solution
csw4 <- eqmcc(ttw4new, details=TRUE, PRI=TRUE, show.cases=FALSE, rowdom=TRUE, use.tilde=FALSE)
csw4

#enhanced parsimonious solution (all simplifying assumptions)

epsw4 <- eqmcc(ttw4new, include="?", details=TRUE, PRI=TRUE, show.cases=FALSE, 
               rowdom=TRUE, use.tilde=FALSE)
epsw4


#Z-Test for enhanced parsimonious solution, 0.8 (almost always sufficient)
# Calculate membership in solution term  tp*op*CM
epsolw4 <- compute("tp*op*CM ", data=tff)

obsepsolw4 <- as.numeric(epsolw4 <= nw)
zepsolw4 <- round(((((sum(obsepsolw4 == 1))/1004) - 0.8)-(1/(2*1004))) / (sqrt((0.8*(1-0.8))/1004)), digits=1)
zepsolw4
pzepsolw4 <- 2*pnorm(-abs(zepsolw4))
pzepsolw4

#intermediate solution with Directional expectations: sp->w, tp->w, op->w.
isw4 <- eqmcc(ttw4new, details = TRUE, include = "?",  use.tilde = FALSE, PRI = TRUE,
              row.dom = TRUE, show.cases = FALSE, 
              dir.exp = "0, 0, 0, -, -")

isw4

#Test for skewness of sufficient conditions (enhanced parsimonious solution)
# calculate membership in paths of M1: tp*op*CM
path1w4 <- compute("tp*op*CM ", data=tff)

A1path1w4 <- round(sum(as.numeric((path1w4 <= nw) & (path1w4 >= tff$W)))/(1004/100), digits = 1)

A2path1w4 <- round(sum(as.numeric((path1w4 >= nw)&(path1w4 >= tff$W)))/(1004/100), digits = 1)

A3path1w4 <- round(sum(as.numeric((path1w4 >= nw)&(path1w4 <= tff$W)))/(1004/100), digits=1)

A4path1w4 <- round(sum(as.numeric((path1w4 <= nw)&(path1w4 <= tff$W)))/(1004/100), digits = 1)


par(mfrow=c(2, 3))
plot(nw, path1w4, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 1',
     xlab= 'Path 1',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path1w4)
text(x=0.8, y=0.5, labels = A2path1w4)
text(x=0.5, y=0.2, labels = A3path1w4)
text(x=0.2, y=0.5, labels = A4path1w4)

A1epsolw4 <- round(sum(as.numeric((epsolw4 <= nw) & (epsolw4 >= tff$W)))/(1004/100), digits = 1)

A2epsolw4 <- round(sum(as.numeric((epsolw4 >= nw)&(epsolw4 >= tff$W)))/(1004/100), digits = 1)

A3epsolw4 <- round(sum(as.numeric((epsolw4 >= nw)&(epsolw4 <= tff$W)))/(1004/100), digits=1)

A4epsolw4 <- round(sum(as.numeric((epsolw4 <= nw)&(epsolw4 <= tff$W)))/(1004/100), digits = 1)


plot(nw, epsolw4, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness eps',
     xlab= 'solution term',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1epsolw4)
text(x=0.8, y=0.5, labels = A2epsolw4)
text(x=0.5, y=0.2, labels = A3epsolw4)
text(x=0.2, y=0.5, labels = A4epsolw4)

#####w5  frequency threshold 10, raw consistency threshold 0.934 #####


ttw5 <- truthTable(tff, outcome="~W", 
                   conditions = "SP, TP, OP, SM, CM", 
                   incl.cut=0.934, n.cut=10, sort.by="incl, n", decreasing=TRUE, 
                   complete=FALSE, show.cases=FALSE)
ttw5

#parsimonious solution
psw5 <- eqmcc(ttw5, include="?", details=TRUE, show.cases=FALSE,  row.dom=TRUE, all.sol=FALSE)

psw5

#exclude contradictory assumptions & rows: OP*CM + SP*TP*CM + sp*TP*OP*SM + SP*TP*op*SM (contradict statement of sufficiency)

ttw5new <- ttw5
suf
rows <- findRows("OP*CM + SP*TP*CM + sp*TP*OP*SM + SP*TP*op*SM", ttw5new, remainders = FALSE)
rows

ttw5new$tt[as.character(rows), "OUT"] <- 0
ttw5new


#conservative solution
csw5 <- eqmcc(ttw5new, details=TRUE, PRI=TRUE, show.cases=FALSE, rowdom=TRUE, use.tilde=FALSE)
csw5

#enhanced parsimonious solution (all simplifying assumptions)

epsw5 <- eqmcc(ttw5new, include="?", details=TRUE, PRI=TRUE, show.cases=FALSE, 
               rowdom=TRUE, use.tilde=FALSE)
epsw5



#Z-Test for enhanced parsimonious solution, 0.8 (almost always sufficient)
# Calculate membership in solution term  sp*TP*op + tp*op*CM + tp*op*SM 
epsolw5 <- compute("sp*TP*op + tp*op*CM + tp*op*SM ", data=tff)

obsepsolw5 <- as.numeric(epsolw5 <= nw)
zepsolw5 <- round(((((sum(obsepsolw5 == 1))/1004) - 0.8)-(1/(2*1004))) / (sqrt((0.8*(1-0.8))/1004)), digits=1)
zepsolw5
pzepsolw5 <- 2*pnorm(-abs(zepsolw5))
pzepsolw5

#intermediate solution with Directional expectations: sp->w, tp->w, op->w.
isw5 <- eqmcc(ttw5new, details = TRUE, include = "?",  use.tilde = FALSE, PRI = TRUE,
              row.dom = TRUE, show.cases = FALSE, 
              dir.exp = "0, 0, 0, -, -")

isw5


#Test for skewness of sufficient conditions (enhanced parsimonious solution)
# calculate membership in paths of M1: sp*TP*op + tp*op*CM + tp*op*SM
path1w5 <- compute("sp*TP*op ", data=tff)
path2w5 <- compute("tp*op*CM ", data=tff)
path3w5 <- compute("tp*op*SM ", data=tff)

A1path1w5 <- round(sum(as.numeric((path1w5 <= nw) & (path1w5 >= tff$W)))/(1004/100), digits = 1)

A2path1w5 <- round(sum(as.numeric((path1w5 >= nw)&(path1w5 >= tff$W)))/(1004/100), digits = 1)

A3path1w5 <- round(sum(as.numeric((path1w5 >= nw)&(path1w5 <= tff$W)))/(1004/100), digits=1)

A4path1w5 <- round(sum(as.numeric((path1w5 <= nw)&(path1w5 <= tff$W)))/(1004/100), digits = 1)


par(mfrow=c(2, 3))
plot(nw, path1w5, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 1',
     xlab= 'Path 1',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path1w5)
text(x=0.8, y=0.5, labels = A2path1w5)
text(x=0.5, y=0.2, labels = A3path1w5)
text(x=0.2, y=0.5, labels = A4path1w5)

A1path2w5 <- round(sum(as.numeric((path2w5 <= nw) & (path2w5 >= tff$W)))/(1004/100), digits = 1)

A2path2w5 <- round(sum(as.numeric((path2w5 >= nw)&(path2w5 >= tff$W)))/(1004/100), digits = 1)

A3path2w5 <- round(sum(as.numeric((path2w5 >= nw)&(path2w5 <= tff$W)))/(1004/100), digits=1)

A4path2w5 <- round(sum(as.numeric((path2w5 <= nw)&(path2w5 <= tff$W)))/(1004/100), digits = 1)

plot(nw, path2w5, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 2',
     xlab= 'Path 2',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path2w5)
text(x=0.8, y=0.5, labels = A2path2w5)
text(x=0.5, y=0.2, labels = A3path2w5)
text(x=0.2, y=0.5, labels = A4path2w5)

A1path3w5 <- round(sum(as.numeric((path3w5 <= nw) & (path3w5 >= tff$W)))/(1004/100), digits = 1)

A2path3w5 <- round(sum(as.numeric((path3w5 >= nw)&(path3w5 >= tff$W)))/(1004/100), digits = 1)

A3path3w5 <- round(sum(as.numeric((path3w5 >= nw)&(path3w5 <= tff$W)))/(1004/100), digits=1)

A4path3w5 <- round(sum(as.numeric((path3w5 <= nw)&(path3w5 <= tff$W)))/(1004/100), digits = 1)

plot(nw, path3w5, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 3',
     xlab= 'Path 3',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path3w5)
text(x=0.8, y=0.5, labels = A2path3w5)
text(x=0.5, y=0.2, labels = A3path3w5)
text(x=0.2, y=0.5, labels = A4path3w5)


A1epsolw5 <- round(sum(as.numeric((epsolw5 <= nw) & (epsolw5 >= tff$W)))/(1004/100), digits = 1)

A2epsolw5 <- round(sum(as.numeric((epsolw5 >= nw)&(epsolw5 >= tff$W)))/(1004/100), digits = 1)

A3epsolw5 <- round(sum(as.numeric((epsolw5 >= nw)&(epsolw5 <= tff$W)))/(1004/100), digits=1)

A4epsolw5 <- round(sum(as.numeric((epsolw5 <= nw)&(epsolw5 <= tff$W)))/(1004/100), digits = 1)


plot(nw, epsolw5, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness eps',
     xlab= 'solution term',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1epsolw5)
text(x=0.8, y=0.5, labels = A2epsolw5)
text(x=0.5, y=0.2, labels = A3epsolw5)
text(x=0.2, y=0.5, labels = A4epsolw5)


#####w6  frequency threshold 10, raw consistency threshold 0.901#####


ttw6 <- truthTable(tff, outcome="~W", 
                   conditions = "SP, TP, OP, SM, CM", 
                   incl.cut=0.901, n.cut=10, sort.by="incl, n", decreasing=TRUE, 
                   complete=FALSE, show.cases=FALSE)
ttw6

#parsimonious solution
psw6 <- eqmcc(ttw6, include="?", details=TRUE, show.cases=FALSE,  row.dom=TRUE, all.sol=FALSE)

psw6

#exclude contradictory assumptions & rows: OP*CM + SP*TP*CM + sp*TP*OP*SM + SP*TP*op*SM (contradict statement of sufficiency)

ttw6new <- ttw6
suf
rows <- findRows("OP*CM + SP*TP*CM + sp*TP*OP*SM + SP*TP*op*SM", ttw6new, remainders = FALSE)
rows

ttw6new$tt[as.character(rows), "OUT"] <- 0
ttw6new


#conservative solution
csw6 <- eqmcc(ttw6new, details=TRUE, PRI=TRUE, show.cases=FALSE, rowdom=TRUE, use.tilde=FALSE)
csw6

#enhanced parsimonious solution (all simplifying assumptions)

epsw6 <- eqmcc(ttw6new, include="?", details=TRUE, PRI=TRUE, show.cases=FALSE, 
               rowdom=TRUE, use.tilde=FALSE)
epsw6



#Z-Test for enhanced parsimonious solution, 0.8 (almost always sufficient)
# Calculate membership in solution term  sp*op + tp*op + op*sm*cm
epsolw6 <- compute("sp*op + tp*op + op*sm*cm ", data=tff)

obsepsolw6 <- as.numeric(epsolw6 <= nw)
zepsolw6 <- round(((((sum(obsepsolw6 == 1))/1004) - 0.8)-(1/(2*1004))) / (sqrt((0.8*(1-0.8))/1004)), digits=1)
zepsolw6
pzepsolw6 <- 2*pnorm(-abs(zepsolw6))
pzepsolw6


#intermediate solution with Directional expectations: sp->w, tp->w, op->w.
isw6 <- eqmcc(ttw6new, details = TRUE, include = "?",  use.tilde = FALSE, PRI = TRUE,
              row.dom = TRUE, show.cases = FALSE, 
              dir.exp = "0, 0, 0, -, -")

isw6

# identify simplifying assumptions and easy counterfactuals

SAw6 <- psw6$SA
SAw6

ECw6 <- isw6$i.sol$C1P1$EC
ECw6

#Z-Test for intermediate solution, 0.8 (almost always sufficient)
# Calculate membership in solution term  sp*op + tp*op + op*sm*cm
isolw6 <- compute("sp*op + tp*op + op*sm*cm ", data=tff)

obsisolw6 <- as.numeric(isolw6 <= nw)
zisolw6 <- round(((((sum(obsisolw6 == 1))/1004) - 0.8)-(1/(2*1004))) / (sqrt((0.8*(1-0.8))/1004)), digits=1)
zisolw6
pzisolw6 <- 2*pnorm(-abs(zisolw6))
pzisolw6

# for path 1
# Calculate membership in path  sp*op 
path1isolw6 <- compute("sp*op ", data=tff)

obspath1isolw6 <- as.numeric(path1isolw6 <= nw)
zpath1isolw6 <- round(((((sum(obspath1isolw6 == 1))/1004) - 0.8)-(1/(2*1004))) / (sqrt((0.8*(1-0.8))/1004)), digits=1)
zpath1isolw6
pzpath1isolw6 <- 2*pnorm(-abs(zpath1isolw6))
pzpath1isolw6

#for path 2
# Calculate membership in path  tp*op
path2isolw6 <- compute("tp*op ", data=tff)

obspath2isolw6 <- as.numeric(path2isolw6 <= nw)
zpath2isolw6 <- round(((((sum(obspath2isolw6 == 1))/1004) - 0.8)-(1/(2*1004))) / (sqrt((0.8*(1-0.8))/1004)), digits=1)
zpath2isolw6
pzpath2isolw6 <- 2*pnorm(-abs(zpath2isolw6))
pzpath2isolw6

# for path 3
# Calculate membership in path  op*sm*cm
path3isolw6 <- compute("op*sm*cm ", data=tff)

obspath3isolw6 <- as.numeric(path3isolw6 <= nw)
zpath3isolw6 <- round(((((sum(obspath3isolw6 == 1))/1004) - 0.8)-(1/(2*1004))) / (sqrt((0.8*(1-0.8))/1004)), digits=1)
zpath3isolw6
pzpath3isolw6 <- 2*pnorm(-abs(zpath3isolw6))
pzpath3isolw6
isw6

#Test for skewness of sufficient conditions (enhanced parsimonious solution)
# calculate membership in paths of M1: sp*op + tp*op + op*sm*cm
path1w6 <- compute("sp*op ", data=tff)
path2w6 <- compute("tp*op ", data=tff)
path3w6 <- compute("op*sm*cm ", data=tff)

A1path1w6 <- round(sum(as.numeric((path1w6 <= nw) & (path1w6 >= tff$W)))/(1004/100), digits = 1)

A2path1w6 <- round(sum(as.numeric((path1w6 >= nw)&(path1w6 >= tff$W)))/(1004/100), digits = 1)

A3path1w6 <- round(sum(as.numeric((path1w6 >= nw)&(path1w6 <= tff$W)))/(1004/100), digits=1)

A4path1w6 <- round(sum(as.numeric((path1w6 <= nw)&(path1w6 <= tff$W)))/(1004/100), digits = 1)


par(mfrow=c(2, 3))
plot(nw, path1w6, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 1',
     xlab= 'Path 1',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path1w6)
text(x=0.8, y=0.5, labels = A2path1w6)
text(x=0.5, y=0.2, labels = A3path1w6)
text(x=0.2, y=0.5, labels = A4path1w6)

A1path2w6 <- round(sum(as.numeric((path2w6 <= nw) & (path2w6 >= tff$W)))/(1004/100), digits = 1)

A2path2w6 <- round(sum(as.numeric((path2w6 >= nw)&(path2w6 >= tff$W)))/(1004/100), digits = 1)

A3path2w6 <- round(sum(as.numeric((path2w6 >= nw)&(path2w6 <= tff$W)))/(1004/100), digits=1)

A4path2w6 <- round(sum(as.numeric((path2w6 <= nw)&(path2w6 <= tff$W)))/(1004/100), digits = 1)

plot(nw, path2w6, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 2',
     xlab= 'Path 2',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path2w6)
text(x=0.8, y=0.5, labels = A2path2w6)
text(x=0.5, y=0.2, labels = A3path2w6)
text(x=0.2, y=0.5, labels = A4path2w6)

A1path3w6 <- round(sum(as.numeric((path3w6 <= nw) & (path3w6 >= tff$W)))/(1004/100), digits = 1)

A2path3w6 <- round(sum(as.numeric((path3w6 >= nw)&(path3w6 >= tff$W)))/(1004/100), digits = 1)

A3path3w6 <- round(sum(as.numeric((path3w6 >= nw)&(path3w6 <= tff$W)))/(1004/100), digits=1)

A4path3w6 <- round(sum(as.numeric((path3w6 <= nw)&(path3w6 <= tff$W)))/(1004/100), digits = 1)

plot(nw, path3w6, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 3',
     xlab= 'Path 3',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path3w6)
text(x=0.8, y=0.5, labels = A2path3w6)
text(x=0.5, y=0.2, labels = A3path3w6)
text(x=0.2, y=0.5, labels = A4path3w6)


A1epsolw6 <- round(sum(as.numeric((epsolw6 <= nw) & (epsolw6 >= tff$W)))/(1004/100), digits = 1)

A2epsolw6 <- round(sum(as.numeric((epsolw6 >= nw)&(epsolw6 >= tff$W)))/(1004/100), digits = 1)

A3epsolw6 <- round(sum(as.numeric((epsolw6 >= nw)&(epsolw6 <= tff$W)))/(1004/100), digits=1)

A4epsolw6 <- round(sum(as.numeric((epsolw6 <= nw)&(epsolw6 <= tff$W)))/(1004/100), digits = 1)


plot(nw, epsolw6, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness eps',
     xlab= 'solution term',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1epsolw6)
text(x=0.8, y=0.5, labels = A2epsolw6)
text(x=0.5, y=0.2, labels = A3epsolw6)
text(x=0.2, y=0.5, labels = A4epsolw6)

#####w7  frequency threshold 50, raw consistency threshold 0.812#####


ttw7 <- truthTable(tff, outcome="~W", 
                   conditions = "SP, TP, OP, SM, CM", 
                   incl.cut=0.812, n.cut=50, sort.by="incl, n", decreasing=TRUE, 
                   complete=FALSE, show.cases=FALSE)
ttw7

#parsimonious solution
psw7 <- eqmcc(ttw7, include="?", details=TRUE, show.cases=FALSE,  row.dom=TRUE, all.sol=FALSE)

psw7

#exclude contradictory assumptions & rows: OP*CM + SP*TP*CM + sp*TP*OP*SM + SP*TP*op*SM (contradict statement of sufficiency)

ttw7new <- ttw7
suf
rows <- findRows("OP*CM + SP*TP*CM + sp*TP*OP*SM + SP*TP*op*SM", ttw7new, remainders = FALSE)
rows

ttw7new$tt[as.character(rows), "OUT"] <- 0
ttw7new


#conservative solution
csw7 <- eqmcc(ttw7new, details=TRUE, PRI=TRUE, show.cases=FALSE, rowdom=TRUE, use.tilde=FALSE)
csw7

#enhanced parsimonious solution (all simplifying assumptions)

epsw7 <- eqmcc(ttw7new, include="?", details=TRUE, PRI=TRUE, show.cases=FALSE, 
               rowdom=TRUE, use.tilde=FALSE)
epsw7


#Z-Test for enhanced parsimonious solution, 0.8 (almost always sufficient)
# Calculate membership in solution term  tp*cm + sp*sm*cm
epsolw7 <- compute("tp*cm + sp*sm*cm", data=tff)

obsepsolw7 <- as.numeric(epsolw7 <= nw)
zepsolw7 <- round(((((sum(obsepsolw7 == 1))/1004) - 0.8)-(1/(2*1004))) / (sqrt((0.8*(1-0.8))/1004)), digits=1)
zepsolw7
pzepsolw7 <- 2*pnorm(-abs(zepsolw7))
pzepsolw7

#intermediate solution with Directional expectations: sp->w, tp->w, op->w.
isw7 <- eqmcc(ttw7new, details = TRUE, include = "?",  use.tilde = FALSE, PRI = TRUE,
              row.dom = TRUE, show.cases = FALSE, 
              dir.exp = "0, 0, 0, -, -")

isw7


#Test for skewness of sufficient conditions (enhanced parsimonious solution)
# calculate membership in paths of M1: tp*cm + sp*sm*cm
path1w7 <- compute("tp*cm ", data=tff)
path2w7 <- compute("sp*sm*cm ", data=tff)

A1path1w7 <- round(sum(as.numeric((path1w7 <= nw) & (path1w7 >= tff$W)))/(1004/100), digits = 1)

A2path1w7 <- round(sum(as.numeric((path1w7 >= nw)&(path1w7 >= tff$W)))/(1004/100), digits = 1)

A3path1w7 <- round(sum(as.numeric((path1w7 >= nw)&(path1w7 <= tff$W)))/(1004/100), digits=1)

A4path1w7 <- round(sum(as.numeric((path1w7 <= nw)&(path1w7 <= tff$W)))/(1004/100), digits = 1)


par(mfrow=c(2, 3))
plot(nw, path1w7, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 1',
     xlab= 'Path 1',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path1w7)
text(x=0.8, y=0.5, labels = A2path1w7)
text(x=0.5, y=0.2, labels = A3path1w7)
text(x=0.2, y=0.5, labels = A4path1w7)

A1path2w7 <- round(sum(as.numeric((path2w7 <= nw) & (path2w7 >= tff$W)))/(1004/100), digits = 1)

A2path2w7 <- round(sum(as.numeric((path2w7 >= nw)&(path2w7 >= tff$W)))/(1004/100), digits = 1)

A3path2w7 <- round(sum(as.numeric((path2w7 >= nw)&(path2w7 <= tff$W)))/(1004/100), digits=1)

A4path2w7 <- round(sum(as.numeric((path2w7 <= nw)&(path2w7 <= tff$W)))/(1004/100), digits = 1)

plot(nw, path2w7, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 2',
     xlab= 'Path 2',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path2w7)
text(x=0.8, y=0.5, labels = A2path2w7)
text(x=0.5, y=0.2, labels = A3path2w7)
text(x=0.2, y=0.5, labels = A4path2w7)

A1epsolw7 <- round(sum(as.numeric((epsolw7 <= nw) & (epsolw7 >= tff$W)))/(1004/100), digits = 1)

A2epsolw7 <- round(sum(as.numeric((epsolw7 >= nw)&(epsolw7 >= tff$W)))/(1004/100), digits = 1)

A3epsolw7 <- round(sum(as.numeric((epsolw7 >= nw)&(epsolw7 <= tff$W)))/(1004/100), digits=1)

A4epsolw7 <- round(sum(as.numeric((epsolw7 <= nw)&(epsolw7 <= tff$W)))/(1004/100), digits = 1)


plot(nw, epsolw7, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness eps',
     xlab= 'solution term',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1epsolw7)
text(x=0.8, y=0.5, labels = A2epsolw7)
text(x=0.5, y=0.2, labels = A3epsolw7)
text(x=0.2, y=0.5, labels = A4epsolw7)

#####w8  frequency threshold 50, raw consistency threshold 0.794 #####


ttw8 <- truthTable(tff, outcome="~W", 
                   conditions = "SP, TP, OP, SM, CM", 
                   incl.cut=0.794, n.cut=50, sort.by="incl, n", decreasing=TRUE, 
                   complete=FALSE, show.cases=FALSE)
ttw8

#parsimonious solution
psw8 <- eqmcc(ttw8, include="?", details=TRUE, show.cases=FALSE,  row.dom=TRUE, all.sol=FALSE)

psw8

#exclude contradictory assumptions & rows: OP*CM + SP*TP*CM + sp*TP*OP*SM + SP*TP*op*SM (contradict statement of sufficiency)

ttw8new <- ttw8
suf
rows <- findRows("OP*CM + SP*TP*CM + sp*TP*OP*SM + SP*TP*op*SM", ttw8new, remainders = FALSE)
rows

ttw8new$tt[as.character(rows), "OUT"] <- 0
ttw8new


#conservative solution
csw8 <- eqmcc(ttw8new, details=TRUE, PRI=TRUE, show.cases=FALSE, rowdom=TRUE, use.tilde=FALSE)
csw8

#enhanced parsimonious solution (all simplifying assumptions)

epsw8 <- eqmcc(ttw8new, include="?", details=TRUE, PRI=TRUE, show.cases=FALSE, 
               rowdom=TRUE, use.tilde=FALSE)
epsw8


#Z-Test for enhanced parsimonious solution, 0.8 (almost always sufficient)
# Calculate membership in solution term  sm*cm
epsolw8 <- compute("sm*cm ", data=tff)

obsepsolw8 <- as.numeric(epsolw8 <= nw)
zepsolw8 <- round(((((sum(obsepsolw8 == 1))/1004) - 0.8)-(1/(2*1004))) / (sqrt((0.8*(1-0.8))/1004)), digits=1)
zepsolw8
pzepsolw8 <- 2*pnorm(-abs(zepsolw8))
pzepsolw8

#intermediate solution with Directional expectations: sp->w, tp->w, op->w.
isw8 <- eqmcc(ttw8new, details = TRUE, include = "?",  use.tilde = FALSE, PRI = TRUE,
              row.dom = TRUE, show.cases = FALSE, 
              dir.exp = "0, 0, 0, -, -")

isw8


#Test for skewness of sufficient conditions (enhanced parsimonious solution)
# calculate membership in paths of M1: sm*cm
path1w8 <- compute("sm*cm ", data=tff)

A1path1w8 <- round(sum(as.numeric((path1w8 <= nw) & (path1w8 >= tff$W)))/(1004/100), digits = 1)

A2path1w8 <- round(sum(as.numeric((path1w8 >= nw)&(path1w8 >= tff$W)))/(1004/100), digits = 1)

A3path1w8 <- round(sum(as.numeric((path1w8 >= nw)&(path1w8 <= tff$W)))/(1004/100), digits=1)

A4path1w8 <- round(sum(as.numeric((path1w8 <= nw)&(path1w8 <= tff$W)))/(1004/100), digits = 1)


par(mfrow=c(2, 3))
plot(nw, path1w8, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 1',
     xlab= 'Path 1',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path1w8)
text(x=0.8, y=0.5, labels = A2path1w8)
text(x=0.5, y=0.2, labels = A3path1w8)
text(x=0.2, y=0.5, labels = A4path1w8)

A1epsolw8 <- round(sum(as.numeric((epsolw8 <= nw) & (epsolw8 >= tff$W)))/(1004/100), digits = 1)

A2epsolw8 <- round(sum(as.numeric((epsolw8 >= nw)&(epsolw8 >= tff$W)))/(1004/100), digits = 1)

A3epsolw8 <- round(sum(as.numeric((epsolw8 >= nw)&(epsolw8 <= tff$W)))/(1004/100), digits=1)

A4epsolw8 <- round(sum(as.numeric((epsolw8 <= nw)&(epsolw8 <= tff$W)))/(1004/100), digits = 1)


plot(nw, epsolw8, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness eps',
     xlab= 'solution term',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1epsolw8)
text(x=0.8, y=0.5, labels = A2epsolw8)
text(x=0.5, y=0.2, labels = A3epsolw8)
text(x=0.2, y=0.5, labels = A4epsolw8)





###############################THEORY EVALUATION COVERAGE ############

tff <- read.csv("tff_d1_c2.csv", row.names=1, header = TRUE, sep = ";",  dec = ",")
describe(tff)
nw <- 1-tff$W



## Theory evaluation for hypothesis 2####

## calculate intersections

Tw <- "sp*tp*op"
Tw
tw <- deMorgan("sp*tp*op")
tw
Sw <- "sp*op + tp*op + op*sm*cm"
Sw
sw <- deMorgan("sp*op + tp*op + op*sm*cm")
sw

TSw <- intersection("sp*tp*op", "sp*op + tp*op + op*sm*cm", snames = "SP, TP, OP, SM, CM")
TSw
factorize (TSw, snames = "SP, TP, OP, SM, CM")

tw
tSw <- intersection("OP + SP + TP", "sp*op + tp*op + op*sm*cm", snames = "SP, TP, OP, SM, CM")
tSw
factorize (tSw, snames = "SP, TP, OP, SM, CM")

sw
tsw <- intersection("OP + SP + TP", "OP + CM*SP*TP + SM*SP*TP", snames = "SP, TP, OP, SM, CM")
tsw
factorize (tsw, snames = "SP, TP, OP, SM, CM")

Tw
Tsw <- intersection("sp*tp*op", "OP + CM*SP*TP + SM*SP*TP", snames = "SP, TP, OP, SM, CM")
Tsw

#calculate membership in intersections

TSw
TSw1 <- compute("sp*tp*op", data=tff)
TSw2 <- compute("sp*tp*op*sm*cm", data=tff)
TSw <- compute("sp*tp*op + sp*tp*op*sm*cm", data=tff)

Tsw

tSw
tSw1 <- compute("SP*tp*op", data=tff)
tSw2 <- compute("SP*op*sm*cm", data=tff)
tSw3 <- compute("sp*TP*op", data=tff)
tSw4 <- compute("TP*op*sm*cm", data=tff)
tSw <- compute("SP*tp*op + SP*op*sm*cm + sp*TP*op + TP*op*sm*cm", data=tff)

tsw
tsw1 <- compute("OP", data=tff)
tsw2 <- compute("SP*TP*OP*CM", data=tff)
tsw3 <- compute("SP*TP*OP*SM", data=tff)
tsw4 <- compute("SP*OP", data=tff)
tsw5 <- compute("SP*TP*CM", data=tff)
tsw6 <- compute("SP*TP*SM", data=tff)
tsw7 <- compute("TP*OP", data=tff)
tsw <- compute("OP + SP*TP*OP*CM + SP*TP*OP*SM + SP*OP + SP*TP*CM + SP*TP*SM + TP*OP", data=tff)



#check for logical remainders
TSw
sum(as.numeric(TSw1 > 0.5))
sum(as.numeric(TSw2 > 0.5))
sum(as.numeric(TSw > 0.5))

sum(as.numeric(tSw1 > 0.5))
sum(as.numeric(tSw2 > 0.5))
sum(as.numeric(tSw3 > 0.5))
sum(as.numeric(tSw4 > 0.5))
sum(as.numeric(tSw > 0.5))

sum(as.numeric(tsw1 > 0.5))
sum(as.numeric(tsw2 > 0.5))
sum(as.numeric(tsw3 > 0.5))
sum(as.numeric(tsw4 > 0.5))
sum(as.numeric(tsw5 > 0.5))
sum(as.numeric(tsw6 > 0.5))
sum(as.numeric(tsw7 > 0.5))
sum(as.numeric(tsw > 0.5))



#calculate % with W and w

TSwW <- round(sum(as.numeric((TSw > 0.5) & (tff$W > 0.5)))/(1004/100), digits = 1)
TSwW

TSww <- round(sum(as.numeric((TSw > 0.5 )&(tff$W < 0.5)))/(1004/100), digits = 1)
TSww

TswW <- round(sum(as.numeric((Tsw > 0.5) & (tff$W > 0.5)))/(1004/100), digits = 1)
TswW

Tsww <- round(sum(as.numeric((Tsw > 0.5 )&(tff$W < 0.5)))/(1004/100), digits = 1)
Tsww

tswW <- round(sum(as.numeric((tsw > 0.5) & (tff$W > 0.5)))/(1004/100), digits = 1)
tswW

tsww <- round(sum(as.numeric((tsw > 0.5 )&(tff$W < 0.5)))/(1004/100), digits = 1)
tsww

tSwW <- round(sum(as.numeric((tSw > 0.5) & (tff$W > 0.5)))/(1004/100), digits = 1)
tSwW

tSww <- round(sum(as.numeric((tSw > 0.5 )&(tff$W < 0.5)))/(1004/100), digits = 1)
tSww



#check for mistakes
fulltew <- sum(TSwW, TSww, tSwW, tSww, TswW, Tsww, tswW, tsww)
fulltew



####Theory evaluation for hypothesis 3#####
TW <- " OP*SM + OP*CM + SP*SM + SP*CM + TP*SM + TP*CM "
TW
tW <- deMorgan("OP*SM + OP*CM + SP*SM + SP*CM + TP*SM + TP*CM ")
tW
SW <- "OP*CM + SP*TP*CM + sp*TP*OP*SM + SP*TP*op*SM"
SW
sW <- deMorgan("OP*CM + SP*TP*CM + sp*TP*OP*SM + SP*TP*op*SM")
sW

TSW <- intersection("OP*SM + OP*CM + SP*SM + SP*CM + TP*SM + TP*CM ", " OP*CM + SP*TP*CM + sp*TP*OP*SM + SP*TP*op*SM ", snames = "SP, TP, OP, SM, CM")
TSW


tW
SW
tSW <- intersection("op*sp*tp + cm*sm", "OP*CM + SP*TP*CM + sp*TP*OP*SM + SP*TP*op*SM", snames = "SP, TP, OP, SM, CM")
tSW

sW
tW
tsW <- intersection("op*sp*tp + cm*sm", "cm*sm + cm*OP*SP + cm*tp + op*sp + op*tp", snames = "SP, TP, OP, SM, CM")
tsW

TW
sW
TsW <- intersection("OP*SM + OP*CM + SP*SM + SP*CM + TP*SM + TP*CM", "cm*sm + cm*OP*SP + cm*tp + op*sp + op*tp", snames = "SP, TP, OP, SM, CM")
TsW


#calculate membership in intersections

TSW
TSW1 <- compute("OP*SM*CM", data=tff)
TSW2 <- compute("SP*TP*OP*SM*CM", data=tff)
TSW3 <- compute("sp*TP*OP*SM", data=tff)
TSW4 <- compute("OP*CM", data=tff)
TSW5 <- compute("SP*TP*OP*CM", data=tff)
TSW6 <- compute("sp*TP*OP*SM*CM", data=tff)
TSW7 <- compute("SP*OP*SM*CM", data=tff)
TSW8 <- compute("SP*TP*SM*CM", data=tff)
TSW9 <- compute("SP*TP*op*SM", data=tff)
TSW10 <- compute("SP*OP*CM", data=tff)
TSW11 <- compute("SP*TP*CM", data=tff)
TSW12 <- compute("SP*TP*op*SM*CM", data=tff)
TSW13 <- compute("TP*OP*SM*CM", data=tff)
TSW14 <- compute("TP*OP*CM", data=tff)

TSW <- compute("OP*SM*CM + SP*TP*OP*SM*CM + sp*TP*OP*SM + OP*CM + SP*TP*OP*CM + sp*TP*OP*SM*CM + SP*OP*SM*CM + SP*TP*SM*CM + SP*TP*op*SM + SP*OP*CM + SP*TP*CM + SP*TP*op*SM*CM + TP*OP*SM*CM + TP*OP*CM", data=tff)

TsW
TsW1 <- compute("SP*OP*SM*cm", data=tff)
TsW2 <- compute("tp*OP*SM*cm", data=tff)
TsW3 <- compute("SP*tp*SM*cm", data=tff)
TsW4 <- compute("SP*tp*op*SM", data=tff)
TsW5 <- compute("SP*tp*op*CM", data=tff)
TsW6 <- compute("SP*TP*OP*SM*cm", data=tff)
TsW7 <- compute("sp*TP*op*SM", data=tff)
TsW8 <- compute("sp*TP*op*CM", data=tff)

TsW <- compute("SP*OP*SM*cm + tp*OP*SM*cm + SP*tp*SM*cm + SP*tp*op*SM + SP*tp*op*CM + SP*TP*OP*SM*cm + sp*TP*op*SM + sp*TP*op*CM", data=tff)

tSW

tsW
tsW1 <- compute("sp*tp*op*sm*cm", data=tff)
tsW2 <- compute("sp*tp*op*cm", data=tff)
tsW3 <- compute("sp*tp*op", data=tff)
tsW4 <- compute("sm*cm", data=tff)
tsW5 <- compute("SP*OP*sm*cm", data=tff)
tsW6 <- compute("tp*sm*cm", data=tff)
tsW7 <- compute("sp*op*sm*cm", data=tff)
tsW8 <- compute("tp*op*sm*cm", data=tff)
tsW <- compute("sp*tp*op*sm*cm + sp*tp*op*cm + sp*tp*op + sm*cm + SP*OP*sm*cm + tp*sm*cm + sp*op*sm*cm + tp*op*sm*cm", data=tff)

#check for logical remainders
TSW
sum(as.numeric(TSW1 > 0.5))
sum(as.numeric(TSW2 > 0.5))
sum(as.numeric(TSW3 > 0.5))
sum(as.numeric(TSW4 > 0.5))
sum(as.numeric(TSW5 > 0.5))
sum(as.numeric(TSW6 > 0.5))
sum(as.numeric(TSW7 > 0.5))
sum(as.numeric(TSW8 > 0.5))
sum(as.numeric(TSW9 > 0.5))
sum(as.numeric(TSW10 > 0.5))
sum(as.numeric(TSW11 > 0.5))
sum(as.numeric(TSW12 > 0.5))
sum(as.numeric(TSW13 > 0.5))
sum(as.numeric(TSW14 > 0.5))

sum(as.numeric(TSW > 0.5))

TsW
sum(as.numeric(TsW1 > 0.5))
sum(as.numeric(TsW2 > 0.5))
sum(as.numeric(TsW3 > 0.5))
sum(as.numeric(TsW4 > 0.5))
sum(as.numeric(TsW5 > 0.5))
sum(as.numeric(TsW6 > 0.5))
sum(as.numeric(TsW7 > 0.5))
sum(as.numeric(TsW8 > 0.5))

sum(as.numeric(TsW > 0.5))

tsW
sum(as.numeric(tsW1 > 0.5))
sum(as.numeric(tsW2 > 0.5))
sum(as.numeric(tsW3 > 0.5))
sum(as.numeric(tsW4 > 0.5))
sum(as.numeric(tsW5 > 0.5))
sum(as.numeric(tsW6 > 0.5))
sum(as.numeric(tsW7 > 0.5))
sum(as.numeric(tsW8 > 0.5))

sum(as.numeric(tsW > 0.5))

#calculate % with W and w

TSWW <- round(sum(as.numeric((TSW > 0.5) & (tff$W > 0.5)))/(1004/100), digits = 1)
TSWW

TSWw <- round(sum(as.numeric((TSW > 0.5 )&(tff$W < 0.5)))/(1004/100), digits = 1)
TSWw

TsWW <- round(sum(as.numeric((TsW > 0.5) & (tff$W > 0.5)))/(1004/100), digits = 1)
TsWW

TsWw <- round(sum(as.numeric((TsW > 0.5 )&(tff$W < 0.5)))/(1004/100), digits = 1)
TsWw

tsWW <- round(sum(as.numeric((tsW > 0.5) & (tff$W > 0.5)))/(1004/100), digits = 1)
tsWW

tsWw <- round(sum(as.numeric((tsW > 0.5 )&(tff$W < 0.5)))/(1004/100), digits = 1)
tsWw

tSWW <- round(sum(as.numeric((tSW > 0.5) & (tff$W > 0.5)))/(1004/100), digits = 1)
tSWW

tSWw <- round(sum(as.numeric((tSW > 0.5 )&(tff$W < 0.5)))/(1004/100), digits = 1)
tSWw



#check for mistakes
fulltew <- sum(TSWW, TSWw, tSWW, tSWw, TsWW, TsWw, tsWW, tsWw)
fulltew


